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Particle  Simulation  of  Auroral  Double  Layers 
Bruce  L.  Smith 

Externally  driven  magnetic  reconnection  has  been  proposed  as  a  possible 
mechanism  for  production  of  auroral  electrons  during  magnetic  substorms.  Fluid 
simulations  of  magnetic  reconnection  lead  to  strong  plasma  flows  towards  the 
increasing  magnetic  field  of  the  earth.  These  plasma  flows  must  generate  large 
scale  potential  drops  to  preserve  global  charge  neutrality.  We  have  examined 
currentless  injection  of  plasma  along  a  dipole  magnetic  field  into  a  bounded 
region  using  both  analytic  techniques  and  particle  simulation. 

Our  analysis  shows  that  the  maximum  potential  for  cold  ions  and  electrons, 
mass  ratio  is  qAV  =  K — oc  /iAB,  where  n  is  their  common  magnetic 
moment,  K  is  the  kinetic  energy  of  the  injected  ions,  and  A B  is  the  difference 
between  the  maximum  magnetic  field  strength,  Bmax  .  and  that  at  injection, 
Bo.  With  thermal  spread  in  particle  magnetic  moments  the  potential  depends  on 
their  temperature  ratios,  and  with  thermal  injection  velocities  the  magnetic 
field  mirror  ratio,  B/gxx  ■  For  drifting  isotropic  Maxwellians,  the  leading  order 

potential  is  qA<j>  =  K~^f ^  g~~~  • 

This  potential  is  a  boundary  condition  for  the  local  field  solution.  While 
global  charge  neutrality  must  be  maintained  in  steady-state,  the  local  poten¬ 
tial  is  determined  by  the  Vlasov-Poisson  equations.  Several  authors  obtained 
quasineutral  solutions  (n*  =  n„)  to  this  boundary  value  problem.  However, 
single-valued  quasineutral  solutions  are  only  accessible  to  special  combinations 
of  fields  and  particle  distributions.  In  our  case  multivalued  solutions  to  the 
quasineutral  equation  arise  for  particular  values  of  electric  and  magnetic  fields. 
A  double  layer  is  that  solution  which  preserves  gross  quasineutrality  within  its 
volume  while  permitting  momentum  balance  between  incident  accelerated  par¬ 
ticles.  The  potential  of  the  double  layer  is  limited  by  the  ion  kinetic  energy  but 
need  not  match  the  global  potential  required  for  overall  charge  neutrality. 

One  and  two  dimensional  electrostatic  particle  simulations  verify  both  dou¬ 
ble  layer  solutions  and  dependence  of  potentials  on  the  injected  energy.  The 
difference  between  global  and  local  potentials  is  absorbed  in  a  sheath  opposite 
the  injection  boundary. 

Potential  formations  are  strictly  dependent  on  the  self-consistent  charge  dis¬ 
tributions  they  support.  Microinstabilities  cause  changes  in  particle  distribu¬ 
tions.  Double  layer  motion  is  associated  with  exchange  of  momentum  between 
particles  and  these  fields. 
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simulations  of  magnetic  reconnection  lead  to  strong  plasma  flows  towards  the 
increasing  magnetic  field  of  the  earth.  These  plasma  flows  must  generate  large 
scale  potential  drops  to  preserve  global  charge  neutrality.  We  have  examined 
currentless  injection  of  plasma  along  a  dipole  magnetic  field  into  a  bounded 
region  using  both  analytic  techniques  and  particle  simulation. 

Our  analysis  shows  that  the  maximum  potential  for  cold  ions  and  electrons, 
mass  ratio  jj,  is  qAV  =  K^—^-  a  fiAB,  where  /z  is  their  common  magnetic 
moment,  K  is  the  kinetic  energy  of  the  injected  ions,  and  A B  is  the  difference 
between  the  maximum  magnetic  field  strength,  Bmax ,  and  that  at  injection, 
B0.  With  thermal  spread  in  particle  magnetic  moments  the  potential  depends  on 
their  temperature  ratios,  p,  and  with  thermal  injection  velocities  the  magnetic 
field  mirror  ratio,  --max. .  For  drifting  isotropic  Maxwellians,  the  leading  order 

potential  is  qA<f>  =  K 1 1^rT*  • 

This  potential  is  a  boundary  condition  for  the  local  field  solution.  While 
global  charge  neutrality  must  be  maintained  in  steady-state,  the  local  poten¬ 
tial  is  determined  by  the  Vlasov- Poisson  equations.  Several  authors  obtained 
quasineutral  solutions  (n,-  =  ne)  to  this  boundary  value  problem.  However, 
single-valued  quasineutral  solutions  are  only  accessible  to  special  combinations 
of  fields  and  particle  distributions.  In  our  case  multivalued  solutions  to  the 
quasineutral  equation  arise  for  particular  values  of  electric  and  magnetic  fields. 
A  double  layer  is  that  solution  which  preserves  gross  quasineutrality  within  its 
volume  while  permitting  momentum  balance  between  incident  accelerated  par¬ 
ticles.  The  potential  of  the  double  layer  is  limited  by  the  ion  kinetic  energy  but 
need  not  match  the  global  potential  required  for  overall  charge  neutrality. 

One  and  two  dimensional  electrostatic  particle  simulations  verify  both  dou¬ 
ble  layer  solutions  and  dependence  of  potentials  on  the  injected  energy.  The 
difference  between  global  and  local  potentials  is  absorbed  in  a  sheath  opposite 
the  injection  boundary. 

Potential  formations  are  strictly  dependent  on  the  self-consistent  charge  dis¬ 
tributions  they  support.  Microinstabilities  cause  changes  in  particle  distribu¬ 
tions.  Double  layer  motion  is  associated  with  exchange  of  momentum  between 
particles  and  these  fields. 
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Externally  driven  magnetic  reconnection  has  been  proposed  as  a  possible 
mechanism  for  production  of  auroral  electrons  during  magnetic  substorms. 
Fluid  simulations  of  magnetic  reconnection  lead  to  strong  plasma  flows 
towards  the  increasing  magnetic  field  of  the  earth.  These  plasma  flows  must 
generate  large  scale  potential  drops  to  preserve  global  charge  neutrality. 
We  have  examined  currentless  injection  of  plasma  along  a  dipole  magnetic 
field  into  a  bounded  region  using  both  analytic  techniques  and  particle 
simulation. 

Our  analysis  shows  that  the  maximum  potential  for  cold  ions  and  elec- 
trons,  mass  ratio  tt,  is  qAV  =  K—^-  <x  (iAB,  where  /i  is  their  common 
magnetic  moment,  K  is  the  kinetic  energy  of  the  injected  ions,  and  A B  is  the 
difference  between  the  maximum  magnetic  field  strength,  Bmax ,  and  that 
at  injection,  Bo.  With  thermal  spread  in  particle  magnetic  moments  the 
potential  depends  on  their  temperature  ratios,  and  with  thermal  injec¬ 
tion  velocities  the  magnetic  field  mirror  ratio,  ■  For  drifting  isotropic 
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Maxwellians,  the  leading  order  potential  is  qA<f>  =  K  b^ax  ' 

This  potential  is  a  boundary  condition  for  the  local  field  solution.  While 
global  charge  neutrality  must  be  maintained  in  steady-state,  the  local  po¬ 
tential  is  determined  by  the  Vlasov-Poisson  equations.  Several  authors 
obtained  quasineutral  solutions  (n,  =  ne)  to  this  boundary  value  problem. 
However,  single- valued  quasineutral  solutions  are  only  accessible  to  special 
combinations  of  fields  and  particle  distributions.  In  our  case  multivalued 
solutions  to  the  quasineutral  equation  arise  for  particular  values  of  electric 
and  magnetic  fields.  A  double  layer  is  that  solution  which  preserves  gross 
quasineutrality  within  its  volume  while  permitting  momentum  balance  be¬ 
tween  incident  accelerated  particles.  The  potential  of  the  double  layer  is 
limited  by  the  ion  kinetic  energy  but  need  not  match  the  global  potential 
required  for  overall  charge  neutrality. 

One  and  two  dimensional  electrostatic  particle  simulations  verify  both 
double  layer  solutions  and  dependence  of  potentials  on  the  injected  energy. 
The  difference  between  global  and  local  potentials  is  absorbed  in  a  sheath 
opposite  the  injection  boundary. 

Potential  formations  are  strictly  dependent  on  the  self-consistent  charge 
distributions  they  support.  Microinstabilities  cause  changes  in  particle  dis¬ 
tributions.  Double  layer  motion  is  associated  with  exchange  of  momentum 
between  particles  and  these  fields. 
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Philosophy 

The  birth  of  plasma  physics  as  a  discipline  is  most  often  associated  with  the 
advent  of  the  fusion  program  in  the  late  50’s.  However,  an  equal  impetus 
for  a  comprehensive  understanding  of  plasma  physics  phenomena  was  gen¬ 
erated  by  coincident  advances  in  space  physics.  Both  the  space  physics  and 
the  fusion  communities  owe  a  great  deal  to  an  even  more  recent  discipline. 

Up  to  a  comparatively  short  time  ago  (1960’s)  the  only  approach  to 
the  many-body  problems  of  plasma  phenomena  was  the  statistical  one  of 
solving  the  coupled  Maxwell-Boltzmann  equations  or  its  moments  -  fluid 
theory  -  applied  to  an  appropriate  model.  Alternatively,  the  experimentalist 
could  construct  an  apparatus  which  closely  matched  the  parameters  under 
investigation  and  with  suitable  diagnostics,  make  desired  measurements.  In 
both  instances  the  construction  of  a  suitable  model,  to  which  the  theory 
or  an  experiment  could  be  tractably  applied,  was  quite  constrained,  and 
certain  simplifying  assumptions  were  necessary  in  order  to  glean  any  useful 
information  at  all. 

With  analytic  theory,  once  approximations  were  made  and  a  solution 
found,  many  of  the  details  of  plasma  behaviour  could  be  obscured.  For 
instance,  reliance  on  Fluid  Theory  provided  no  information  about  wave- 
particle  interactions.  Similarly,  experiments  were  limited,  in  that  not  all 
parameters  could  be  measured.  These  measurements  depended  severely  on 
the  ingenuity  and  pursestrings  of  the  experimenter.  In  many  respects  then 
the  most  useful  results  of  theory  and  experiment  are  to  demonstrate  very 
general  phenomena  while  application  to  specific  phenomena  were  limited 
to  embellishing  the  general  results  with  appropriate  detail  -  as  for  example, 
explanation  of  Whistler  waves  in  the  ionosphere. 


ix 


With  advent  of  powerful  computers  a  radical  new  approach  to  plasma 
physics  became  possible.  Instead  of  statistical  solution  to  the  many-body 
problem,  one  could  actually  follow  the  self-consistent  trajectories  of  each 
particle  in  an  ensemble.  The  data  from  each  particle  could  be  stored  and 
used  to  provide  an  “exact”  description  of  particles  and  fields. 

Of  course,  this  is  an  oversimplification.  Computer  time  and  memory  are 
limited  and  other  complications  abound.  First,  the  particles  can  only  be 
followed  in  discrete  steps  so  that  a  series  of  snapshots  and  not  a  continuous 
motion  picture  is  possible.  Second,  calculations  of  the  fields  in  a  pair-wise 
manner  was  prohibitively  time  expensive  and  memory  consuming.  Third, 
the  number  of  particles  which  could  be  followed,  and  therefore  particle 
density,  was  generally  much  lower  than  that  of  the  modeled  phenomena. 
Finally,  since  the  model  is  necessarily  finite,  questions  about  what  to  do 
with  particles  and  fields  at  the  boundaries  is  a  thorny  one. 

All  of  these  considerations  mean  that,  while  a  powerful  new  tool  ex¬ 
ists  to  examine  specific  physical  situations,  the  role  of  the  modeler,  the 
computational  physicist,  is  an  enormous  one.  Once  the  model  has  been 
constructed  and  sucessfully  implemented  for  a  particular  situation,  it  is 
important  to  know  that  the  model  can  be  “scaled”  to  a  broader  range  of 
additional  problems. 

Even  more  important  than  the  code  as  implemented  is  the  ability  to 
organize  and  use  a  set  of  “tools”  into  working  code.  So,  while  useful  results 
may  be  obtained  using  another’s  code,  its  construction  or  even  its  modifi¬ 
cation,  clearly  relies  on  the  experience  gained  by  doing.  This  experience  is 
ultimately  the  most  important  product  of  this  effort. 
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Chapter  1 
Introduction 


Since  1958  when  the  Van  Allen  belts  were  discovered,  we  have  come  to  rec¬ 
ognize  complex  interactions  between  the  sun  and  the  earth  coupled  by  the 
sun’s  plasma  (the  solar  wind),  its  radiation  (at  the  earth’s  ionosphere)  and 
their  magnetic  fields  (the  earth’s,  the  Interplanetary  Magnetic  Field  (IMF), 
and  current-driven  magnetic  disturbances.)  One  of  the  most  intriguing  and 
beautiful  manifestations  of  these  interactions  is  the  aurora  borealis. 

As  pictured  in  the  accompanying  sketch  (Figure  1.1),  observation  of 
these  aurorae  is  not  new.  Greeks  recorded  them  many  centuries  ago  but 
their  name  is  attributed  to  Gassandi  (1621)[12,  p.15].  The  detailed  physical 
processes  which  lead  to  aurorae,  however,  are  still  being  understood.  The 
cause  of  the  emissions  themselves  weren’t  known  until  the  50 ’s  when,  from 
electron  densities  at  the  aurorae,  they  were  linked  to  impact  of  energetic 
electrons  upon  oxygen  molecules  in  the  earth’s  ionosphere. 

Since  their  discovery  two  different  types  of  aurorae  have  been  identi¬ 
fied.  (Figure  1.2).  Diffuse  aurorae  are,  just  as  their  name  implies,  broad 
bands  of  light  which  appear  almost  continuously  in  the  polar  region  sky. 
The  broad  distribution  of  electron  energies  which  cause  diffuse  aurorae  are 
associated  with  pitch  angle  scattering  into  the  electron  loss  cone  above  the 
earth’s  dipole  magnetic  field. 

In  contrast,  discrete  (or  Bright  Active)  aurorae,  are  identified  as  ribbons 
of  light.  In  1960  Mcllwain  [31,  p.99]  discovered  a  characteristic  electron 
energy  peak  at  6  keV  during  discrete  aurorae.  For  these  energetic  electrons 
to  reach  the  earth,  accelerating  potentials  must  exist  parallel  to  the  earth’s 
magnetic  field  [31,  p.  99ffj  and  a  source  of  energy  must  be  available  to  drive 
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Figure  1.1:  A  woodcut  by  Fridtjof  Nansen;  Nansen  depicts  himself  strolling 
on  the  ice  under  a  triple  curtain-like  form  of  the  aurora;  the  auroral  arcs. 
(From  Nansen’s  Nord  I  Takeheunen,  1911)  [2,  Cover  Page]. 
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Figure  1.2:  Schematic  diagram  showing  the  main  characteristics  of  auroral 
displays,  as  seen  from  above  the  north  geomagnetic  pole.  Auroral  arcs  are 
indicated  by  lines  and  the  diffuse  aurora  is  indicated  by  the  shaded  region. 
This  is  slightly  modified  from  the  original  version  given  by  Akasofu  (1976) 
[2.  P-4]. 

them. 

Discrete  aurorae  are  correlated  with  increased  solar  activity  and  the 
direction  of  the  IMF  with  respect  to  the  earth’s  dipole  magnetic  field.  These 
in  combination  are  believed  to  cause  substorms.  During  these  substorms 
plasma  is  observed  to  flow  from  the  tailwaxd  region  toward  and  away  from 
the  earth.  As  modeled  by  Sato  [38]  the  energy  involved  in  this  flow  results 
from  externally  driven  magnetic  reconnection. 

A  “semi-empirical”  theory  to  substantiate  plasma  flows  as  a  source  of 
energetic  electrons  for  substorm-generated  aurorae  is  due  to  Serizawa  and 
Sato[42].  They  developed  an  expression  for  the  potential  drop  between  in¬ 
jection  and  loss  boundaries  required  to  preserve  global  charge  neutrality. 
Their  analysis  revealed  that  for  ion  parallel  kinetic  energy,  K,  and  temper¬ 
ature  ratio,  T{/Te ,  e<f> max  ~  AT/(1  -f  Ti/Te )  with  only  mild  dependence  on 
electron  to  ion  mass  ratios  m/M  <€.  1  and  mirror  ratios  B max /Bo  1. 
As  noted  by  the  authors,  this  theory  yields  no  details  about  the  nature  of 
background  particle  interaction  with  the  local  potential. 

Under  the  assumption  of  quasineutrality  several  authors  have  shown, 
both  in  theory  [34], [52]  and  through  numerical  calculation  [14],  that  large, 
long- scale  length  ( l  Ad)  potential  drops  along  magnetic  field  lines  may 
develop  self-consistently.  Although  this  mechanism  has  been  proposed  as 
a  candidate  for  acceleration  of  auroral  electrons,  results  between  differ- 
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ent  models  are  contradictory  and  heavily  dependent  on  assumed  particle 
distributions.  Further,  the  assumption  of  quasineutrality  limits  acquired 
solutions  to  only  a  subset  of  those  actually  possible. 

Among  alternatives  proposed  for  electron  acceleration  are  double  layers 
(DLs)  -  potential  drops  over  relatively  short  scale  lengths  ~  ^).  Both 
strong  ( e<{>  kT )  and  weak  (e<f>  ~  kT)  DLs  have  recently  been  observed 
above  the  earth’s  poles  [49]. 

The  role  of  double  layers  in  laboratory  and  space  plasmas  has  been  the 
subject  of  much  investigation.  In  an  early  analysis  Block  [6,  p.349]  modelled 
four  species  of  particles-reflecting  and  passing,  electrons  and  ions-incident 
upon  a  strong  DL.  Two  fluid  equations,  an  adiabatic  equation  of  state,  and 
Poisson’s  equation  lead  to  criteria  on  the  drift  velocities  of  ions  and  electrons 
incident  on  the  high  and  low  potential  sides  of  the  DL,  respectively.  These 
are  called  Bohm  criteria  in  analogy  to  similar  criteria  on  ions  in  a  plasma 
sheath[7,  p.77].  An  analgous  Langmuir’s  condition  [30,  p.  973]  leads  to  the 
perceived  requirement  for  net  current  in  all  DLs. 

Focussing  on  these  criteria  as  a  recipe,  one  could  easily  construct  a 
DL.  Plasma  experiments  and  simulations  required  only  fixed  potentials  at 
the  boundaries  to  drive  required  currents  or  a  floating  potential  (or  even 
periodic  boundary  conditions)  with  large  enough  drifts  between  species. 
Because  of  these  drifts,  these  plasmas  were  inherently  unstable.  (See  Fig. 
1.3.) 

Just  recently  double  layers  have  been  demonstrated  under  more  relaxed 
assumptions  about  boundary  conditions  and  particle  distributions.  Some 
authors  have  found  ways  to  relax  the  constraints  imposed  by  the  Bohm  cri¬ 
teria.  For  example,  Kan  and  Lee  [25]  concluded  that  the  condition  on  the 
electron  velocity  is  obviated  by  presence  of  trapped  electrons.  In  their  sup¬ 
port  simulations  by  Wagner  showed  double-layers  to  develop  in  the  presence 
of  just  a  current  sheet[51]. 

Sato  and  Okuda  [36,37]  performed  a  series  of  simulations  using  even 
“more  realistic”  conditions.  Their  model  (Fig.  1.4(a))  was  that  of  polar 
region  field  lines  in  a  self-consistent  circuit  (Fig.  1.4(b)).  Initial  conditions 
included  a  driving  potential  and  current.  The  subsequent  potential  and 
current  were  related  by  a  fixed  resistance  consistent  with  the  initial  con¬ 
ditions.  Specifically,  they  assumed  an  electron  drift  velocity,  v*  <  vthe, 
and  Te  »  T<.  This  is  the  range  of  parameters  for  ion- acoustic  instabili¬ 
ties,  but  avoids  the  large  relative  drifts  which  may  cause  other  two-stream 
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Figure  1.3:  Initial  propagation  of  a  pulse  driven  by  fixed  potential  boundary 
conditions.  Top  left:  the  potential  profile.  Top  right:  electron  distribution 
function  at  the  high-potential  boundary.  Bottom  left:  electron  phase  space. 
Bottom  right:  ion  phase  space.  (After  Hubbard  and  Joyce  [35].)  [43] 
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Figure  1.4:  Model  of  polar  region  current  driven  by  the  convection  elec¬ 
tric  field  (a).  The  plasma  is  embedded  in  a  self-consistent  circuit  (b)  [37, 
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Figure  1.5:  Sato  and  Okuda  results  exhibit  multiple  weak  DLs.  These  DLs 
move  with  a  velocity  near  the  ion-acoustic  velocity  but  recur  so  that  the 
number  of  DLs  is  approximately  constant[37,  p.3364]. 

instabilities. 

One  of  the  results,  shown  in  Figure  1.5,  was  obtained  for  Vdjvthe  —  0-6. 
As  apparent,  the  simulation  resulted  in  multiple  weak  ( ecf>  ~  Te)  DLs  about 
1000 Ad  apart  and  with  scale  lengths  1  ~  50 Ad-  These  DLs  are  unstable 
and  propagate  at  near  the  ion-acoustic  velocity  but  recur  at  a  rate  such 
that  the  number  of  DLs  is  approximately  constant. 

Hasegawa  and  Sato  [22]  provided  the  theoretical  mechanism  for  such 
DLs.  Basically,  an  ion  hole  is  created  which  cuts  off  the  electron  current. 
Formation  of  an  adjacent  electron  hole  follows.  This  yields  a  DL  which 
decays  on  the  ion  time  scale. 

Most  recently  laboratory  and  theory  experiments  have  demonstrated 
double  layer  solutions  even  under  currentless  conditions.  Stenzel  et  al  [44] 
conducted  an  experiment  in  which  a  dipole  magnetic  field  reflected  an  in¬ 
cident  ion  beam.  This  experiment  resulted  in  inherently  currentless  strong 
DLs  for  varying  magnetic  field  strengths.  Cohen  et  al  modeled  the  likely 
formation  of  double  layers  in  the  throat  of  a  mirror  device[15].  Perkins  and 
Sun  [35]  generalized  their  results  by  demonstrating  eigenvalue  solutions  for 
currentless  DLs. 
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Interestingly,  quasineutrality  does  not  preclude  formation  of  double  lay¬ 
ers  (DLs)  [52,  p.  1526] [45].  These  results  contrast  with  those  like  Chiu  and 
Schultz’  [14],  using  the  condition  of  local  charge  neutrality  to  demonstrate 
the  requirement  for  particular  mixes  of  particle  velocity  distributions,  but 
agree  with  more  recent  work  by  Stern,  who  includes  double  layers  within 
an  overall  quasi-neutral  framework. 

As  described  in  this  introduction  it  is  our  goal  to  extend  work  on  dou¬ 
ble  layers  to  more  realistic  scenarios  such  as  these  “recent”  theories,  ex¬ 
periments,  and  simulations.  Chapter  2  provides  a  brief  survey  of  auroral 
physics,  leading  to  a  physical  context  for  construction  of  our  model.  We  are 
motivated  by  observations  that  aurorae  are  caused  by  energetic  electrons. 
In  particular,  we  examine  one  possible  source  of  these  electrons  suggested  by 
Sato  et  al  [38],  [42]  that  during  substorms  plasma  flows  earthward  driven 
by  magnetic  reconnection  in  the  tail  region  leading  to  large  scale  poten¬ 
tials  to  be  explored  by  plasma  constituents.  Their  results  are  supported 
by  globed  calculations,  however,  and  rely  on  static  injection  distributions. 
They  therefore  provide  few  details  about  the  local  behaviour  of  particles 
and  the  self-consistent  electric  potential  in  space  and  time.  Alternative  so¬ 
lutions  exist  which  satisfy  the  boundary  conditions  of  the  global  potential 
drop,  injected  particle  distributions  and  background  magnetic  field. 

In  Chapter  3  we  review  theories  and  measurements  from  experiments 
and  observations  relevant  to  auroral  double  layers,  and  other  plasmas  im¬ 
mersed  in  a  dipole  magnetic  field.  Many  of  the  previous  theories  suffer 
from  the  limited  focus  on  global  versus  local  analysis  of  potential  solutions 
or  ignore  entirely  the  possiblity  for  existence  of  double  layers  by  excluding 
the  class  of  solutions  to  which  they  belong.  These  deficiencies  are  partially 
due  to  unwarranted  regard  for  thresholds  and  criteria  rather  than  the  un¬ 
derlying  physics  from  which  they  were  developed.  This  chapter  has  the 
simultaneous  goal  of  exposing  the  reader  to  analytic  and  graphical  tech¬ 
niques  for  analysing  solutions  to  Poisson’s  equation.  Use  of  both  global 
charge  conservation,  requiring  currentless  boundari°s,  and  local  quasineu¬ 
trality,  exploring  solutions  to  Poisson’s  equation  by  mapping  An  =  0,  are 
accomodated  by  the  same  physics. 

We  have  examined  solutions  to  the  Vlasov- Poisson  equations  using  both 
analytic  and  computational  models.  Among  the  questions  we  ask:  What 
is  the  local  distribution  of  potential-double  layer  or  otherwise?  May  they 
be  currentless  or  do  criteria  prohibit  formation  of  double  layers?  Are  they 
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accessible?  Do  they  persist  or  are  they  subject  to  instabilities?  How  do 
global  and  local  approaches  coincide? 

In  Chapter  4  we  obtain  both  the  global  and  local  analytic  solutions. 
These  analyses  weave  together  the  theories  of  the  previous  chapter  while 
developing  observations  and  supporting  the  model  assumptions  for  the  sim¬ 
ulations  to  follow.  While  Serizawa  and  Sato  use  “semi-empirical”  methods 
to  obtain  a  global  requirement  for  large  scale  potential  drops,  we  develop 
an  analytic  result.  We  analyze  expected  simulation  results  for  the  local 
potential  using  quasi-neutral  techniques.  In  particular  we  assume  Drifting 
Maxwellian  particle  distributions,  beginning  in  the  limit  of  cold  temper¬ 
atures  for  both  Tj_  and  Tj|  and  proceeding  to  finite  values  for  both.  The 
cold  distribution  is  used  explicitly  by  Schmidt  [40]  as  an  example  of  mov¬ 
ing  plasma  behaviour,  but  is  implicit  in  the  assumed  distributions  of  other 
authors  [35],  [15].  We  find  that  the  underlying  global  solutions  retain  their 
dependence  on  the  injected  ion  kinetic  energy  even  for  finite  T±  and  7]|  but 
local  DL  solutions  exist  only  for  warm  Tj|  electrons. 

None  of  the  analytic  techniques  permit  observation  of  temporal  and 
spatial  dependence.  Particle  simulations  allow  observation  of  both,  thus 
demonstrating  accessibility,  the  role  of  instabilities,  interaction  betweeen 
the  self-consistent  particle  distributions  and  how  local  requirements  on  the 
potential  satisfy  competing  global  boundary  conditions.  They  permit  ob¬ 
servation  beyond  physically  limited  experiments  (such  as  those  in  space) 
and  give  the  “experimenter”  more  freedom  in  terms  of  flexibility  and  avail¬ 
ability  of  “apparatus”  and  expenditure  of  time  and  money.  Data  can  be 
stored  and  additional  “non-intrusive”  diagnostics  performed.  These  diag¬ 
nostics  are  limited  only  by  the  ingenuity  of  the  computational  physicist  and 
simulation  time  and  memory  available. 

Previous  simulations  were  inadequate  since  they  concentrated  on  peri¬ 
odic,  non-bounded  simulation  regions  or  fixed  potentials  and  lacked  a  real¬ 
istic  magnetic  field  and  injection  scheme.  In  Chapter  5,  using  particle  simu¬ 
lation,  we  modeled  a  flowing  neutral  plasma  injected  along  a  (fully)  dipole 
magnetic  field.  The  computational  model  is  that  of  a  three-dimensional 
magnetic  mirror  with  a  one-dimensional  (thus  ID)  electrostatic  solution 
to  Maxwell’s  equations  lying  along  the  mirror  axis.  A  currentless,  neutral 
(one-sided)  Maxwellian  plasma  is  injected  at  the  low-field  (x=L)  side  (e.g., 
magnetospere)  into  an  initially  vacuum  simulation  region  and  propagates 
toward  the  high-field  (x=0)  side  (e.g.,  the  ionosphere.)  Field  boundary 
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conditions  are  ^>(0)  =  0,  d<f>/dx(L )  =  — E(L )  =  0[4S] .  Model  diagnostics 
include  phase-space  plots;  local  charge  densities,  fields,  and  potentials;  fre¬ 
quency  and  spatial  spectra;  and  energy  calculations.  A  variety  of  particle 
distributions,  boundary  conditions  and  equations  of  motion  were  used. 

Both  double  layers  and  large  scale  potentials  with  eV  kT  are  ob¬ 
served.  The  double  layers  move  in  correspondence  with  changes  in  velocity 
space  distributions  of  incident  particles.  These  changes  coincide  with  ob¬ 
served  instabilities.  Movement  of  the  double  layer  in  time  is  associated 
with  Langmuir’s  condition.  Comparisons  are  made  between  theory  and 
“experiment.” 

Extension  of  the  simulation  to  two- dimensions  offers  the  possibility  to 
observe  instabilities,  transport  processes  and  scale  lengths  not  permitted 
in  ID,  but  modeling  in  2D  requires  additional,  sometimes  subtle,  consid¬ 
erations.  In  chapter  6  we  describe  construction  of  and  results  from  a  two- 
dimensional  code.  Comparison  is  made  to  behaviour  in  ID. 

In  Chapter  7  we  conclude  by  summarizing  our  contributions,  suggesting 
future  employment  of  our  models  and  identify  gaps  yet  to  be  explored. 

General  methods  of  particle  simulation  are  described  in  a  set  of  appen¬ 
dices  to  detail  techniques  and  tools  necessary  for  code  construction.  Impor¬ 
tant  difficulties  had  to  be  overcome  in  the  simulation  of  a  plasma  injected  at 
a  non-periodic  boundary.  In  ID  we  developed  techniques  to  combat  alias¬ 
ing  of  Fourier  transform  solutions.  In  2D  we  extended  use  of  Buneman’s 
algorithm  [11]  for  solving  Poisson’s  equation  to  mixed  boundary  conditions. 
In  both  codes  we  developed  and  used  a  guiding  center  pusher  for  electrons 
which  permitted  use  of  larger  time  steps.  For  this  reason,  although  we 
initiated  construction  of  these  simulations  with  a  particular  application  in 
mind,  we  emphasize  these  general  aspects  to  relate  the  model’s  construction 
to  possible  follow-on  investigation. 

Finally,  an  extensive  bibliography  should  add  to  the  usefulness  of  this 
document  as  a  reference  source. 
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Chapter  2 
Aurorae  1 


2.1  The  Sun’s  Plasma  and  Magnetic  Fields 

The  aurorae  result  from  complex  plasma  phenoma  under  the  influence  of 
interactions  between  the  sun  and  the  earth.  Regardless  of  their  proximity, 
these  plasmas  remained  relatively  inaccessible  until  comparatively  recently 
and  their  study  is  a  relatively  new  branch  of  space  physics.  In  this  chapter 
we  introduce  the  physics  of  these  plasmas  which  provide  motivation  for  our 
dissertation  and  a  context  within  which  to  build  analytic  and  computational 
models. 

2.1.1  The  Sun 

The  Sun  is  the  major  source  of  energy  and  a  major  source  of  particles 
for  interactions  manifest  in  phenomena  such  as  aurorae.  The  core  at  the 
sun’s  center,  heated  by  gravitational  pressure,  produces  over  1026  W  via  p-p 
and  carbon  cycle  fusions.  This  energy  radiates  out  through  the  radiation 
zone  producing  X-rays,  boils  up  through  the  convection  zone  and  is  visibly 
apparent  in  the  photosphere.  In  this  transit  the  temperature  decreases  from 
107oK  at  the  core  to  5  x  103oK  at  the  photosphere.  However,  in  the  outer 
chromosphere  large  amplitude  acoustic  or  shock  waves  heat  the  gas  of  the 

1  Except  where  otherwise  noted,  the  substance  of  this  chapter  was  gleaned  from  the 
1986  Princeton  University  AS556  Lecture  Notes,  Selected  Topics  in  Space  Physics,  by 
Professor  Hideo  Okuda  [33].  Author  references  with  years  in  parentheses  are  from  these 
notes. 
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photosphere  to  106oK  and  densities  continually  decrease  from  1017  cm-3 
primarily  neutrals  to  a  plasma  density  of  106  —  109  cm-3  at  the  corona. 

The  corona  is  the  location  of  many  solar  activities.  Among  these  are  so¬ 
lar  flaxes  and  sunspots.  Sunspots  are  dark  spots  (low  temperature  regions) 
on  the  sun’s  surface  which  are  also  observed  to  be  local  regions  of  strong 
magnetic  fields  (~  lkG)  much  larger  than  the  normal  sun  magnetic  field,  a 
dipole  of  approximate  magnitude  B=1  Gauss  with  a  22  year  reversal  cycle. 
Sunspots  are  thought  fo  originate  from  the  wrapping  of  magnetic  field  lines 
as  the  sun  rotates  (Babcock,  1959).  This  model  also  describes  the  11  year 
sunspot  cycle  with  its  familiar  “butterfly  pattern.” 

Solar  flares  violently  release  magnetic  energy  associated  with  these  sunspots 
as  the  field  lines  buoy  up  above  the  coronal  surface.  The  energy  stored  in 
solar  flares  is  as  much  as  1032  ergs  and  can  be  released  in  tens  of  minutes. 

Its  release  may  be  explained  by  “driven”  (Pettschek-Parker)  magnetic  field 
line  reconnection.  This  mechanism  may  also  be  important  in  the  Earth’s 
magnetic  tail. 

Radiation  from  the  sun  may  be  approximated  as  blackbody  at  5800° K 
and  is  the  primary  portion  of  energy  released  by  the  sun.  At  the  earth  this 
radiation  is  reponsible  for  visible  light-day  and  night-as  well  as  production 
of  the  ionosphere,  where  the  high  frequency  edge  (TJV  and  X-rays)  of  the 
radiated  spectrum  is  absorbed. 

2.1.2  The  Solar  Wind 

The  solar  wind  was  first  postulated  by  Biermann  (1951)  to  explain  the 
anti-sunward  orientation  of  comets  tails.  Radiation  pressure  alone  was 
insufficient  to  account  for  this  phenomena.  Observations  have  since  shown 
that  typical  earth  parameters  for  the  solar  wind  are  n=5  cm-3,  T=2  x  105oK 
and  u=400  km/sec.  The  wind  consists  of  90%  H  and  10%  He  with  a  mean 
free  path  lmfp  ~  1  AU  =  150  x  106  km  (500  light  sec).  The  energy  flux 
from  the  solar  wind  is  1013  —  1014  W  and  is  the  driver  for  magnetospheric 
activity.  One  may  readily  verify  that  the  kinetic  energy  density  of  the  solar 
wind  at  the  earth  is  much  greater  than  its  magnetic  field  energy  density. 

A  theoretical  foundation  for  the  solar  wind  was  advanced  by  Chapman 
(1958).  His  was  a  static  model  in  which  the  pressure  of  the  corona  was 
balanced  by  the  gravitational  force.  This  model  was  unsuccessful  in  that 
the  gravitational  potential  was  too  weak  to  contain  the  corona.  An  alternate 


12 


dynamic  model,  in  which  outgoing  solar  wind  is  replaced  by  solar  gas,  was 
more  successful.  A  particular  manifestation  of  this  model  is  supersonic 
expansion  of  the  solar  wind  ( Ma  =  5)  in  the  gravitational  throat  of  the 
sun. 


2.1.3  Interplanetary  Magnetic  Field 

The  dipole  field  of  the  sun  is  dragged  along  by  the  solar  wind.  Because  of 
the  sun’s  rotation,  the  field  lines  form  an  Archimedes  spiral.  At  the  earth 
the  resulting  magnetic  field  points  at  an  angle  \  ~  45°  to  the  line  between 
the  sun  and  earth  in  the  sun’s  equatorial  plane.  The  magnitude  of  the  field 
at  the  earth  is  ~  5  x  lO-5^?  and  for  distances  much  greater  than  the  earth’s 
is  primarily  in  the  <f>  direction.  At  large  distances  from  the  sun  density 
decreases  as  1/r2  while  Bimf  oc  1/r.  outward  limit  for  the  solar  wind 
is  obtained  where  magnetic  pre  ssure  balances  particle  pressure  at  50  AU, 
well  beyond  Pluto. 


2.2  The  Earth’s  Magnetosphere  (See  Fig.  2.1. 

2.2.1  The  Earth’s  Magnetic  Field 

For  up  to  about  10  earth  radii  the  magnetic  field  of  the  earth  is  well  modeled 
by  a  simple  magnetic  dipole  produced  by  an  infinitesimally  small  current 
loop.  Such  a  field  is  characterized  by  its  magnetic  moment  (m)  and  is  given 
by  the  formula  [24,  p.  182] 


B  = 


3n(n  ■  m)  —  m 


Since  the  earth  m  =  —mz,  B  may  be  well  expressed  in  cylindrical  coordi¬ 
nates 

s  m  3pz  3z! 

b  =  ^I-7T?  +  (1-—  M 

and  its  magnitude  easily  computed  as 

B  =  ^(1  +  ^)1/2 
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Magnatopaust 


If  the  field  is  known  at  a  given  point,  say  x0,  then  in  terms  of  that  field 

The  magnitude  of  the  earth’s  field  at  1  earth  radius  is  known  to  be  .312  G 
at  the  equator  (and  .628  G  at  the  poles.)  The  total  field  B  =  Bb  is  easily 
seen  to  be 


3 pz  ^  3z2 


and  the  two  directions  of  b  =  bpp  +  bzz  are  given  by 

K  =  -3g(l  +  3£)-''J 

i,  =  (1  -  3fJ)(l  +  3JI-)-1/3 


While  the  IMF  represents  but  a  small  perturbation  to  the  near  earth 
magnetic  field.  Bimf  is  approximately  that  of  the  earth  at  r  ~  (^io-i  )1/3  ~ 
15  earth  radii. 


2.2.2  The  Magnetosheath 

The  first  encounter  of  the  solar  wind  with  the  earth’s  magnetic  field  is  at 
the  bow  shock  (R  ~  15  Re)  where  the  supersonic  flow  from  the  sun  de¬ 
celerates  to  subsonic  with  respect  to  va  .  The  solar  wind  plasma  heats  in 
the  magnetosheath  exchanging  kinetic  for  thermal  energy,  and  is  finally 
blocked  at  the  magnetopause  (R  ~  10  Re)  where  the  solar  wind  pressure 
balances  pressures  from  the  earth’s  magnetic  field  and  particles  in  the  mag¬ 
netosphere,  within  which  the  earth’s  magnetic  field  is  confined.  The  distant 
extent  of  the  magnetosheath  is  R  ~  200  Re  where  the  magnetic  fields  of 
the  earth  find  the  solar  wind  become  comparable[31,  p.86]. 

As  the  plasma  flows  past  the  earth,  it  drags  along  the  earth’s  magnetic 
field  creating  two  categories  of  field  lines-open  and  closed.  Subject  to  adi¬ 
abatic  constraints  particles  are  free  to  enter  or  leave  the  earth’s  ionosphere 
primarily  at  the  polar  caps  where  the  magnetic  field  lines  are  open.  In 
absence  of  collisions  or  wave-particle  interactions  the  solar  wind  may  not 
penetrate  the  closed  lines  of  the  magnetosphere.  However,  two  popular 
mechanisms  have  been  proposed  to  allow  for  this  occurence. 
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Figure  2.2:  Mapping  of  the  convection  electric  potentials  and  the  surface 
S  from  the  magnetosphere  to  the  ionosphere,  as  seen  from  above  the  north 
pole  [31,  p.  60]. 

The  first,  proposed  by  Axford-Hines,  invokes  the  Kelvin-Helmholtz  in¬ 
stability  as  a  source  of  “anomalous  viscosity.”  In  this  case  the  solar  wind 
may  penetrate  the  earth’s  magnetic  field  on  a  steady  state  basis.  Addi¬ 
tionally,  in  periods  of  increased  solar  activity,  when  the  IMF  may  acquire 
random  components  in  the  N-S  direction,  Dungey  [31,  p.56-59]  theorized 
that  a  southward  component  would  cause  reconnection  both  sunward  and 
anti-sunward  of  the  earth,  permitting  free  access  of  the  solar  wind  into  the 
earth’s  magnetosphere. 

2.2.3  Convection  Electric  Field 

As  the  solar  wind  traverses  the  earth’s  magnetic  field,  an  electric  field  is 
mapped  into  the  ionosphere  along  open  dipole  magnetic  field  lines.  This 
field  drives  currents  perpendicular  to  the  magnetic  field  and  establishes 
current  patterns  which  map  into  polar  caps.  The  observed  value  of  this 
electric  field  is  20 mV/m  over  the  polar  caps,  amounting  to  an  integrated 
potential  of  e<f>  ~  50  kV[31,  p.64].  In  the  polar  region  ionosphere  this 
potential  causes  the  flow  of  currents  both  perpendicular  to  the  magnetic 
field  and  parallel  (Pederson  currents)  and  perpendicular  (Hall  currents)  to 
the  electric  field[31,  p.60].  Birkeland  currents  flow  parallel  to  the  magnetic 
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field. 


2.2.4  Ionosphere 

Nearest  the  earth  the  ionosphere  covers  the  globe.  Because  of  solar  radi¬ 
ation  and  neutral  density  fall  off  with  height,  different  regions  of  plasma 
may  be  identified.  Below  the  E  region  (~  100  km)  the  collision  frequency 
is  such  that  uen  >  Qe  and  i/,„  >  fij,  above  it  Q,e  >  uen  >  i/tn  >  fi,-.  As  a 
result,  the  E  region  is  the  site  of  current  amplification,  and  Hall  conductiv¬ 
ity  dominates  current  flow  in  the  ionosphere.  Ionospheric  plasma  may  exit 
through  the  polar  caps  creating  the  polar  wind. 

2.2.5  Plasma  Sheet 

Directly  behind  the  earth  the  ambient  polar  wind  collects  to  form  the 
plasma  sheet.  The  plasma  sheet  is  the  demarcation  between  the  two  oppo¬ 
sitely  directed  components  of  the  combined  earth-IMF  magnetic  field.  The 
sheet  is  populated  with  plasma  particles  T  =  .1  —  10fceV,n  =  .01  —  1cm-3 
and  V  =  10  —  1000fcm/s[31,  p.2]. 

Because  of  the  oppositely  directed  magnetic  fields  a  current  must  flow 
in  the  neutral  sheet  region  defined  by  B=0.  This  is  known  as  the  sheet 
current.  The  plasma  sheet  is  a  tremendous  reservoir  of  solar  wind  energy 
E  ~  3  x  1022  —  1025  ergs[31,  p.3j.  Part  of  the  sheet,  known  as  the  tail, 
stretches  more  than  200  Re  as  estimated  by  mapping  polar  electric  fields 
into  the  tail  [31,  p.86]. 

Examination  of  individual  particle  orbits  as  they  enter  the  neutral  sheet 
show  that  the  particles  are  trapped  and  oscillate  between  opposite  mag¬ 
netic  poles  of  the  plasma  sheet  fields.  Electrons  and  ions  drift  in  opposite 
directions-the  resulting  current  is  in  the  direction  of  the  electric  field.  The 
work  done  on  the  particles  is  P~  1012  W[31,  p.87]. 

Because  of  its  current  system,  the  sheet  is  particularly  vulnerable  to 
instabilities  from  reconnection  (or  tearing  modes.)  These  instabilities  make 
the  sheet  a  source  for  other  of  the  earth’s  regions. 
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2.2.6  Plasmasphere 

Just  outside  the  earth  is  the  closed  region  known  as  the  earth’s  plasmas¬ 
phere.  It  is  a  cold  (T  <  1  eV)  [31,  p.2],  dense  (n  ~  103)  region  of  plasma 
which  ends  abruptly  at  the  plasmapause  where  the  corotation  field  of  the 
earth,  Ec  =  —vr  x  B/c,  coupled  to  the  convection  electric  field  sweeps  out 
the  region  between  Re= 3-6  [31,  p.76]. 

Throughout  the  region  between  the  ionosphere  and  the  sheet  are  the 
famous  ring  currents  and  radiation  belts  where  high  energy  electrons  and 
ions  are  trapped.  The  particles  are  responsible  for  a  magnetic  field  depres¬ 
sion,  AJ3  ~  lOOsnT  [31,  p.3],  at  the  earth’s  surface  during  magnetic  storms. 
This  region  also  stores  a  significant  amount  of  energy-2  x  1022— 1025  ergs[31, 
p.  2].  The  rings  are  accelerated  by  ExB,  Fermi  and  betatron  acceleration 
and  are  depopulated  by  ion  cyclotron  waves  and  charge  exchange  [31,  p.  3]. 

2.2.7  Plasma  Parameters 

Parameters  throughout  these  regions  are  conveniently  summarized  in  Fig. 
2.3.  Using  n=l-100  cm"3,  B=102  -  104  x  10"5G,  and  Te  <  T,  ~  100's  eV 
yields  u}pe  ~  uce  ~  105  —  106  rad/s  and  j3  <  1.  In  this  parameter  regime 
the  electrostatic  approximation  is  appropriate  [28]. 


2.3  Aurorae  and  Substorms 

2.3.1  Diffuse  Aurorae 

The  lights  known  as  aurorae  are  caused  by  energetic  electrons  striking  am¬ 
bient  molecules,  primarily  oxygen  at  the  4756A  line,  within  _the  earth’s 
ionosphere.  Within  the  earth’s  magnetic  field  an  electric  field  E  =  —v  x  B 
forms  to  achieve  a  fluid  force  balance.  The  mapping  of  this  field  back  into 
the  tail  drives  the  tail  plasma  earthward. 

Electrostatic  waves  such  as  ES  electron  cyclotron  harmonics  have  been 
shown  to  be  responsible  for  scattering  (primarily)  convected  electrons  into 
the  loss  cone  causing  the  diffuse  aurora.  The  aurorae  display  character¬ 
istic  patterns  associated  with  return  currents  required  to  replenish  these 
electrons. 
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Figure  2.3:  Model  of  the  near- Earth  plasma  in  equatorial  regions  find  at 
middle  latitudes;  R/Rq  is  the  geocentric  distance  in  Earth  radii  [3]. 
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Figure  2.4:  Illustration  of  convective  flow  V  in  the  equatorial  plane  resulting 
from  the  superposition  of  the  interplanetary  and  the  earth’s  magnetic  field. 
The  surface  S,  charged  as  indicated,  separates  the  regions  of  sunward  and 
anti-sunward  convection[31,  p.59]. 

2.3.2  Discrete  Aurorae 

The  visible  forms  of  discrete  aurorae  have  scale  widths  of  ~  10 km  and 
are  imbedded  in  larger  scale  ~  100 Arm  inverted  V  structures  explained  by 
current  patterns  established  above.  However,  the  current  densities  as  mea¬ 
sured  for  discrete  aurorae,  j  ~  10-6  —  10 ~5A/m2,  cannot  be  explained  by 
flows  into  a  normal  loss  cone  [31,  p.  104],  The  loss  cone  must  be  widened  by 
parallel  electric  fields  with  integrated  potential  ~  few  keV.  Possible  mech¬ 
anisms  for  these  fields  include  large  scale  potentials  generated  by  suitable 
particle  distributions  or  both  strong  and  multiple  weak  DLs.  Backscattered 
particles  and  Auroral  Kilometric  Radiation  (AKR)  from  the  RH  extraor¬ 
dinary  wave,  rivaling  Jupiter  as  a  source  of  radio  emissions,  are  associated 
with  these  potential  structures[31,  p.145]. 

2.3.3  Geomagnetic  Substorms 

During  increases  in  solar  activity  changes  in  the  solar  wind  and  IMF  take 
place.  The  plasma  from  the  solar  wind  may  collect  at  the  plasma  sheet 
more  rapidly  than  it  dissipates  due  to  normal  diffusion  processes.  Because 
the  magnetic  field  is  so  low,  the  resulting  system  may  be  unstable,  injecting 
plasma  along  field  lines  toward  the  polar  caps  to  cause  increased  auroral 
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Figure  2.5:  Summary  of  the  distribution  and  directions  of  the  large-scale 
field-aligned  currents  as  determined  from  Triad  satellite  magnetic  field  ob¬ 
servations.  The  hatched  area  near  noon  indicates  confused  current  di¬ 
rections  (from  Iijima  and  Potemra,  1976b)  (©by  American  Geophysical 
Union)[31,  p.105]. 


21 


Upward  Accelerated 
Positive  Ions  ■ — 

(Beams,  Conics) 


Kilometric  Radiation 


Downward  Accelerated 
Electrons  \ 


1/  J 


Electric 

Equipotentials 


Ion  Conics- 


Upward 

Field- 

Aligned 

Currents 


ih 


// 


- VLF  Hiss 

— VLF  Saucers 

Downward  Field- 
Aligned  Currents 


Aurora 


Ionosphere 


■Magnetic  Field  Line 


22 


activity  known  as  auroral  storms.  During  these  auroral  storms,  the  ambient 
densities  may  increase  above  typical  values. 

Substorms  are  a  subset  of  a  storm  and  last  for  several  hours  while  a 
storm  may  last  several  days.  Substorms  are  associated  with  discrete  aurorae 
and  are  evidenced  by  BEN  (Broadband  Electrostatic  Noise)  and  ion  conics. 

The  actual  process  for  formation  of  a  substorm  is  controversial.  Sub¬ 
storms  may  be  the  result  of  local  tearing  mode  instabilities  and  can  cause  a 
heightened  flow  of  plasma  toward  the  earth  as  well  as  a  tailward  flow  known 
as  a  plasmoid.  Parameters  for  this  flow  are  n=.l-l  cm-3,  V=100-1500  km/s, 
and  Te  <  Ti  <  2keV.  A  substorm,  related  to  tearing  mode  instabilities,  is 
thought  to  be  driven  by  some  external  source  and  this  process  is  therefore 
called  “driven  reconnection.” 


2.3.4  Plas^n  Jet  Theory 

Using  a  2D  r '  ID  model  with  anomalous  resistivity,  Sato  simulated  an 
externally  driven  magnetic  reconnection  creating  strong  plasma  jets  with 
speed*  vA  «  500 km/s  [38].  Subsequently  Sato  and  Hasegawa  [39]  de¬ 
veloped  a  theory  for  this  process.  While  tearing  modes  saturate  at  low 
flow  velocity,  externally  driven  flows  have  large  velocities  consistent  with 
observation  of  strong  jetting  of  protons  [38,  p.  7178]. 

A  subsequent  model  of  potentials  generated  from  these  flows  is  proposed 
as  a  mechanism  for  auroral  electron  acceleration.  Serizawa  and  Sato  showed 
that  an  equilibrium  solution  exists  which  is  currentless  and  requires  no 
trapped  particles.  To  do  this,  they  specified  the  particle  distribution  at  an 
injection  point  into  a  mirror  field  and  then  used  the  curentless  condition, 
ji  =  jei  to  obtain  the  required  maximum  potential  at  the  opposite  boundary. 

The  injected  distributions  axe 


/;(t>||0,O  =  Nj(—)*e^~ajV,±o)  x  {e[-Q>Ho-“;)2]  -  el— }  (2.1) 


where  j=e,i  and  a  =  (^),  u  is  the  flow  speed  and  N  the  particle  density.  For 
frequencies  less  than  the  ion  cyclotron  frequency  but  greater  than  collision 
frequencies  we  may  assume  that  both  energy  and  magnetic  moment  are 
conserved.  From  these  assumptions  the  return  flux  at  the  injection  point 
may  be  calculated.  Imposition  of  j,  =  je  then  yields  a  solution  for  the 


23 


Figure  2.7:  Sketches  illustrate  a  model  of  the  magnetotail  structure  change 
associated  with  magnetic  reconnection.  (Top)  The  normal  state.  (Center) 
A  state  of  the  magnetotail  just  after  magnetic  reconnection  sets  in.  In 
this  stage  ,  accelerated  plasmas  on  the  earth  side  of  the  reconnection  point 
are  still  trapped,  and  those  on  the  antisolar  side  are  confined  in  a  mag¬ 
netic  island.  (Bottom)  Illustration  of  a  stage  of  the  substorm  expansion  in 
which  trapped  particles  axe  precipitated  down  into  the  ionosphere  and  the 
magnetic  island  is  propelled  tailward  [38,  p.  7178]. 
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maximum  potential  which  may  be  obtained  numerically  as 


e<f>MAX  =  -A/uq(1  - 


BmaxAT1 


(2.2) 


This  scenario  will  serve  as  a  model  for  our  simulations  and  theory.  Al¬ 
though  this  formalism  does  not  specify  the  particle  densities  or  potential  at 
any  point  between  the  source  point  and  throat,  it  provides  a  useful  formula 
with  which  to  compare  our  results.  We  shall  devote  the  remainder  of  the 
dissertation  to  examining  both  global  and  local  solutions  for  the  potential 
and  considering  the  relationship  between  these  approaches.  Additionally, 
this  technique  for  obtaining  the  global  potential  will  be  used  in  Chapter  4 
to  analyze  requirements  for  global  charge  neutrality. 
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Chapter  3 

Sheaths  and  Double-Layers 


The  previous  chapter  introduced  a  method  for  determining  global  potential 
requirements  for  a  bounded  region.  In  this  chapter  we  consider  techniques 
for  analyzing  local  self-consistent  field/particle  distributions  and  the  result¬ 
ing  observations  from  their  application. 


3.1  Maxwell-Boltzmann  Equations 


Electromagnetic  phenomena  are  governed  by  Maxwell’s  equations: 


V  •  E  =  47 rp 

3.1(a) 

VxE  =  - 

3.1(6) 

V  ■  B  =  0 

3.1(c) 

vxj  =  j(4  *;+#) 

3.1(d) 

These  equations  are  familiar  to  us  as  Gauss’  Law  (3.1(c)),  Poisson’s  Equa¬ 
tion  (3.1(a)),  Faraday’s  Law  (3.1(b))  and  Ampere’s  Law  (3.1(d)).  Gauss’ 
Law  and  Poisson’s  Equation  may  be  treated  as  initial  value  conditions 
which,  once  satisfied,  remain  so  by  virtue  of  Faraday’s  and  Ampere’s  Laws. 
It  is  often  useful  to  define  the  electric  field  in  terms  of  a  scalar  and  vector 
potentials 


IdA 

c  dt 


with  V  •  A  —  0,  so  that  Poisson’s  equation  may  be  equally  well  expressed 
as  —  V2<£  =  4t rp. 
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The  charge  and  current  densities  p  and  j  can  be  obtained  by  summing 
the  zeroth  and  first  velocity  moments  of  the  particle  distribution  function, 
(/),  for  each  plasma  component  (o),  weighted  by  its  charge  (q): 


P  —  ^aQa  J  dvfc  —  EiaqaTla 

j  —  ^ otQoi  j  J a  —  ^Of^faTar 


The  particle  distributions  in  turn  evolve  according  to  the  Boltzmann 
equation, 

+  v  *  V /  +  a  •  Vvf  —  —-^\coiii  (3-2) 

where  the  RHS  of  (3.2)  represents  “collisions”,  essentially  the  particle/field 
interaction  terms  not  included  in  its  LHS.  When  the  term  on  the  RHS  of 
(3.2)  is  small  and  collisions  may  be  ignored,  we  have  the  Vlasov  equation, 

^  +  v-  V/  +  a  ■  V„/  =  0  (3.2a) 

The  acceleration,  a,  is  commonly  given  in  terms  of  the  Lorentz’  force, 

F  =  ma  =  q(E  +  v  x  B)  (3.3) 


and  couples  these  Maxwell-Boltzmann  equations. 


3.2  The  Fluid  Equations 

The  fluid  equations  derive  from  velocity  moments 

J  dviF{---) 

of  Boltzmann’s  Equation.  These  equations  are  not  closed,  since  each  in¬ 
volves  one  higher  moment.  “Closure”  is  obtained  by  truncating  the  series 
at  some  moment  and  replacing  it  with  some  reasonable  form. 

As  presented  above,  one  set  of  equations  may  be  obtained  for  each 
plasma  constituent.  The  fluid  equations  for  a  plasma  of  electrons  and  one 
specie  of  ions  are  commonly  known  as  two-fluid  theory  and  are  familiar  as 

equation  of  continuity:  |y  +  V  •  nv  =  0  (3.4a) 

equation  of  motion:  mn[§7  +  (^  *  V)u]  =  qn(E  +  v  x  B)  —  Vp  (3.46) 
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where  we  have  assumed  a  scalar  pressure,  p.  Often,  the  isotropic  ideal  gas 
law,  p=nkT,  may  be  used.  E  and  B  obtain  from  Maxwell’s  equations  (3.1) 
above. 

Alternatively,  by  combining  (adding  or  subtracting)  these  individual 
species  equations,  one  obtains  one-fluid  theory.  These  are 

conservation  of  mass:  +  V  •  (/>mu)  =  0  (3.5a) 

conservation  of  charge:  §f  +  V  •  j  =  0  (3.5 b) 

equation  of  motion:  pm  §7  =  jxB  —  Vp  +  pmg  (3.5 c) 

Generalized  Ohm’s  Law:  E  +  v  x  B  =  ijj  +  ^(j  x  B  —  Vpe)  (3.5 d) 

In  the  above,  mass  flow  is  dominated  by  the  ions  while  electrical  properties 
tend  to  depend  on  the  electron’s  motion.  The  resistivity  rj  is  defined  by  the 
first  velocity  moment  of  the  collision  term  in  Boltzmann’s  equation  (3.2). 


3.3  Useful  Idealizations 

3.3.1  Quasi-Neutrality  (Plasma  Approximation) 

By  nature  a  plasma  tends  toward  charge  neutrality.  (Opposite  charges 
attract,  like  charges  repel.)  A  useful  approximation  is  to  assume  that  the 
densities  of  ions  and  elections  are  equal 

An  =  n,  —  ne  =  0  (3-6) 

This  assumption  is  known  as  the  plasma  approximation  or  quasineutrality. 
Of  course,  with  this  assumption  the  value  of  the  electric  field  can  no  longer 
be  freely  determined  directly  from  Poisson’s  equation  but  must  instead 
come  from  some  other  means  such  as  the  fluid  equations  (3.4)  and  (3.5). 

It  is  useful  to  examine  the  validity  of  the  plasma  approximation.  We 
do  so  in  one-dimension.  For  a  single  ion  species  Poisson’s  equation  (3.1(a)) 
can  be  rewritten  by  changing  variables  to 

(PV 

=  Nt-N,  (3.7) 

where  A  &  =  -J  4^7,  N  =  n/n0,  n0  =  n,  +  ne,  V  =  —e<f>/kT,  and  X  =  x/Ap. 

(See  [43,  eq(14)].)  Since  |iVe|  and  |iV,|  <  1,  jiV,-  —  ATe|  <  1  and,  averaging 
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Poisson’s  equation  over  a  scale  length,  L, 


1_ 

L 


Ni  -Ne><  1 


where  <  •••  >  denotes  the  average  j;  fo  dX(-  •  •).  Quasineutrality  holds 
when  \A&$\  <  1  or  the  scale  length  for  the  change  in  slope  of  the  poten¬ 
tial  is  much  greater  than  a  Debye  length.  This  is  trivially  true  for  |||r|  C  1 
over  the  same  interval.  Although  quasineutrality  may  not  hold  locally,  it 
must  hold  in  the  macro  sense,  even  in  the  formation  of  sheaths  and  DLs[45]. 


3.3.2  The  Adiabatic  Approximation 

The  steady-state  Vlasov  equation  may  be  easily  solved  for  the  case  of  one¬ 
dimensional  electrostatic  fields  and  a  mirror  force.  When  the  magnetic 
moment,  /x,  is  an  adiabatic  invariant,  the  force  on  the  particles  may  be 
expressed  in  terms  of  a  scalar  potential  (4>  =  q<j>  +  fiB): 


d$  d<f>  dB 
=  =  ~Qd^  +  fl'cb 

Entering  this  for  the  acceleration  term,  the  steady-state  Vlasov  equation 
becomes 

ya/_  j_d$a/  =  o 

dx  m  dx  dv 

and  may  be  rewritten  as 

dV  dK 

where  V  =  $/kT  and  K  =  1/2 mv2/kT. 

Treating  V  and  K  as  coordinates  in  two  space,  we  see  that 


=  0 


V/  •  (i  -  j)  =  0 

where  V  =  Vi,  K  —  Kj  and  V  =  We  may  represent  this  equation 

graphically  as  a  function  (f)  whose  gradient  is  orthogonal  to  the  direction 
of  i  —  j.  (See  Figure  3.1.)  The  solution  to  this  equation  corresponds 
to  rotating  the  axes  45°  such  that  V/  lies  perpendicular  to  the  family 
of  curves  H=K-fV= const  ant.  Therefore  f  is  a  member  of  the  family  of 
functions,  f(H+c),  c  constant.  The  reader  may  recognize  lines  parallel  to 
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Figure  3.1:  K-V  Plane  showing  lines  of  constant  K  ±  V.  These  lines  may 
be  identified  with  the  Hamiltonian,  H=K+V,  and  the  Lagrangian,  L=K-V. 


V/  as  defining  a  set  of  Lagrangians,  L  =  K-V  =  c12i...,  while  H  may  be 
identified  as  the  Hamiltonian. 

If  we  use  the  boundary  condition  that  f  is  a  drifting  Maxwellian  in  the 
absence  of  a  potential  (at  V=0),  we  know  that 

f(K  +  c)  =  /*, 

and  c  just  determines  the  constant  of  normalization.  To  maintain  this  func¬ 
tional  form,  f(H)  =  This  is  simply  the  Boltzmann  or  adiabatic 

distribution.  For  a  cold  species  (T=0)  the  contribution  to  the  Maxwellian 
occurs  only  at  £=constant.  This  is  equivalent  to  conservation  of  energy. 

A  more  familiar  version  of  the  adiabatic  approximation  is  obtained  from 
the  first  law  of  thermodynamics  when  we  insist  that  an  ideal  gas  remains 
isolated  (TdS=0)  during  a  process  (See  Fig.  3.2.) 


TdS  = 
dW  =  pdV  = 

U  = 

—kTd  In  n  = 
n  = 


dW  +  dU  =  0 

p(d—),  p=nkT  (ideal  gas) 
n 

N(]-mu 2  +q<f>  +  pB )  =  £ 
£ 

d{^mu2  +  q<f>  +  pB) 
n0e~e^kT 
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dQ  =  Tds 


Figure  3.2:  Diagram  depicting  flow  of  energy  into  and  out  of  an  isolated 
control  volume. 


In  practice  the  adiabatic  approximation  is  used  when  the  distribution 
has  time  to  reach  thermodynamic  equilibrium-when  we  know  that  we  have  a 
Maxwell-Boltzmann  distribution-but  before  species  have  time  to  exchange 
energy  with  other  species  (or  waves.) 
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3.4  Sheaths 


3.4.1  One  -  Dimensional  Sheath 

The  solution  for  the  electric  potential  in  a  plasma  of  warm  electrons  and 
cold  ions  in  contact  with  a  conducting  wall  is  attributed  to  Langmuir.  This 
geometry  is  depicted  at  Figure  3.3.  To  the  left  of  x  =  0  we  take  E  =  0, 
implying  n,  =  ne  =  no.  We  use  the  adiabatic  equation  of  state  for  electrons, 

ne  =  n0e$, 

and  combine  the  equation  of  continuity  (3.4a) 

V • (nv)  =  0 
(nv)  o  =  (nv) 

and  conservation  of  energy 

^Md2  +  e<f>  =  \Mv  o 

v  = 


to  obtain  the  ion  particle  density 

ni  =  no(1  ~  p£|r’ 

Substituting  into  Poisson’s  equation  (3.1(a))  and  defining  V  =  we 
obtain 

£H„( _ l _ e-) 

where  the  scale  length  and  debye  length  are  defined  above  and  the  Mach 
number  is  defined  as  M  =  in  terms  of  the  sound  speed  c,  =  (^j)?. 
Multiplying  by  =  £  and  using  £0  =  0  leads  to: 

\lx^  “  +  £>**•  + ‘"I 
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Figure  3.3:  The  potential  <f>  in  a  planar  sheath.  Cold  ions  are  assumed  to 
enter  the  sheath  with  a  uniform  velocity  uo.  The  regions  of  the  sheath  are 
marked  as  the  presheath,  sheath  and  anode  (d)  [13,  p.245]. 


This  equation  can  only  be  satisfied  for 


9V  i 

1  +  M2  <  (1  +  +  e~v 

M.1 

when  £2  >  0.  For  small  F,  near  X  =  0  [13,  p.246], 


1 

2 


F2(l 


1 


M2 

so  M2 


)  > 


> 

v0  >  ca 


0 

1  or 


Equivalently,  we  may  require  for  self-consistency  that  at  V=0  >  0  and 

£  >  o, 


2Fr3/2J_ 

M2’  M2 


+exp_K1_ 
exp  — F 
M2 


vo 


> 

> 

> 

> 


0 

(1  + 

1 


2V  s-3/2  1 

M2’  M2 


c. 


as  before.  This  condition  is  known  as  the  Bohm-Sheath  criterion  and  spec¬ 
ifies  a  minimum  ion  injection  velocity,  v0,  for  the  formation  of  sheaths. 
Alternatively,  this  criteria  has  been  shown  above  as  a  condition  on  the 
slope  of  the  relative  charge  density.  We  will  exploit  this  equivalence  in 
observations  to  follow. 

The  demarcation  between  the  sheath  and  the  presheath  is  marked  by  the 
existence  of  a  quasineutral  or  “plasma”  solution  (eqn.  3.6).  This  is  evident 
in  Langmuir’s  evaluation  of  low-pressure,  long  mean  free  path  plasma,  as 
pictured  in  Fig.  3.4.  The  reader  should  note  that  the  formal  plasma  solution 
is  double- valued  in  <f>  and  no  unique  “plasma”  solution  exists  beyond  s  =  s0. 

Outside  the  sheath  region  in  the  presheath,  v  is  axbitraxy.  To  get  a 
velocity  v  >  ct,  there  must  be  an  accelerating  electric  field  in  the  presheath 
plasma  while  still  satisfying  the  condition  of  quasineutrality.  This  condition 
may  still  be  satisfied  provided  the  scale  length  is  long.  The  field  itself  comes 
from  the  distribution  of  ions  and  electrons  produced  within  the  plasma. 
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Figure  3.4:  The  relation  of  the  plasma  and  complete  solutions  of  the 
plasma-sheath  integral  equation  [50,  Fig.  6,  p.902].  (Here  tj  and  s  cor¬ 
respond  to  the  variables  <f>  and  x.) 
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At  some  point  near  the  wall  the  contribution  of  electrons  to  the  space 
charge  becomes  negligible.  Ignoring  this  contribution,  we  calculate  a  po¬ 
tential  as  measured  from  the  sheath  boundary. 

. _ /,  e<t> 


2 dx^dx 


=  -[8™„MU’{[1  -  M]i  -  1}  _ 


For  $$  >  1  and  (§*)’  «  8jrn„MU5(|^)i, 

d<j>  ,  .  2c(f> .  I,  i 

=  [8-n.MU„(-— )>]» 


.aw3'" 


=  [87rn0Mu0(2e/M)2]2 


|$|3/4  =  3/4[87rnoMuo(2e/M)5]?x 
For  a  thickness  d,  the  current  density,  j0,  is  related  to  <f>  by 

3 

jo  =  4/9(2e/M)4^j  (3.8) 

This  is  known  as  the  Child-Langmuir  law  for  space  charge  limited  flow[13, 
p.  248].  jo  is  established  such  that  global  charge  neutrality  is  maintained 
and  is  fixed  by  the  ion  production  rate,  fa  is  such  that  the  flow  of  electrons 
is  equal  to  the  ion  flow  rate,  d  then  changes  to  suit  these  two  parameters.  A 
similar  relationship  is  obtained  for  a  cold,  zero  field  cathode  where  d  is  the 
distance  between  the  cathode  (injection  point)  and  the  anode  at  potential 
<t>d .  We  shall  exploit  these  facts  in  the  analyses  that  follow. 

3.5  Double  Layers  with  Current 

3.5.1  Double  Sheaths 

A  similar  analysis  to  the  last  can  be  made  for  the  case  of  cold  ions  and 
electrons  incident  upon  a  potential  drop  within  a  plasma.  This  analysis 
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0 


1 


Figure  3.5:  A  generic  model  for  double  layer  analysis  with  passing  and 
trapped  ions  and  electrons  in  each  region  0  and  1.  The  height  of  the  lines 
indicate  the  relative  velocities  of  charges  in  response  to  the  potential,  <j>. 

involves  participation  of  four  different  species-  passing  and  trapped,  elec¬ 
trons  and  ions.  (See  Figure  3.5.)  The  following  set  of  equations  relate  the 
particle  boundary  conditions  to  the  potential: 
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Flux: 

Energy: 


(nu)i  =  (nu)io  ,  (nu)e  —  (nu)e, 
=  e(4>o  ~  <t>)  ,  ^  =  e<f> 

-  4>)]'i  ,  «e=(?)^ 

Poisson:  ^  (n^  —  ne)e 


/o^f  =  -fidxfrwe- 


§■  =  \frn/2e[jeV^  +  \jM/mji(y/<f)o  -  <j>  -  y^o)] 

If  E  =  0  at  <£  =  0O>  ie  =  This  is  Langmuir’s  condition  for  a  double 

sheath.  Noting  from  the  last  section  je  =  \p£ji  for  single  sheaths,  this 
analysis  is  equivalent  to  treating  two  separate  plasma  sheaths  in  contact. 

In  this  form  Langmuir’s  condition  clearly  represents  a  pressure  balance 
between  the  electric  field  on  the  LHS  and  the  particle  impulse  on  the  other. 
Because  the  momentum  gained  by  the  ions  and  electrons  in  the  potential 
drop  4> o  are  proportional  to  the  square  root  of  their  mass  ratios,  the  rate  at 
which  they  can  be  provided  to  the  double  sheath,  their  current  densities, 
are  inversely  so.  In  general,  Langmuir’s  condition  is  given  by  the  expression 
[30,  Eqn.  (49),  p.  973] 

r<t>  o 

/  pd<t>  =  0  (3.9) 

when  the  electric  field  at  either  end  point  is  negligible. 


3.5.2  Block  Theory 

Block  undertook  a  more  complete  study  of  DLs.  He  derived  conditions  for 
ionospheric  DLs  under  the  condition  of  frozen-in  magnetic  flux.  Ignoring 
collisions  and  gravity,  he  derived  self-consistency  requirements, 

.r  _  ...  dB 

(mu2)i0  >  ')TiQ  +  Te0{—— — ^§-)*=0 


and  (mu2)el  >  jTel  +  7h(  — — i-Ly§-)*=i> 


eE  +  MeO-57 
,eE  +  ph* , 

eE  —  fin 


where,  here,  7  is  the  ratio  of  specific  heats.  (In  terms  of  degrees  of  freedom 
d,  7  =  2±i.)  These  two  conditions  are  Bohm  criteria.  We  observe  that,  for 


dB 

dz 


<  0  and  E  >  0,  the  criteria  can  be  satisfied  for  any  velocities  provided 


and 


e|g|+W>lgl  >  Ik 
e|£|-/.«lfl  “  7T« o 

e|£|  +  /*«.lfl  Ik 
elEI-W.lfl  “  7ra 


However,  when  mu{ £  <  yTi,e  +Te,i,  these  constraints  limit  the  electric  field. 
Therefore,  these  stronger  criteria  hold  for  strong  DLs. 

These  Bohm  criteria  provided  strict  conditions  for  DL  formation-that 
is,  the  directed  velocities  of  the  ions  and  electrons  on  each  side  of  the  DL 
must  be  relatively  great.  Since  the  current  density  is 


j  =  e[(nu),0  +  (mi)ei], 


these  criteria  led  to  the  common  observation  that  currentless  double  layers 
could  not  exist. 

Similarly,  Block  developed  a  Langmuir  condition  identical  to  that  of  the 
double  sheath  case  in  the  last  section,  =  (jg-)*.  However,  Block  hinted 
at  another  possibility  for  the  existence  of  double  layers, 


The  Langmuir  condition  requires  a  supply  of  electrons  and 
ions  in  the  right  proportion  from  both  sides  of  the  double  layer. 

If  that  is  not  possible,  the  layer  will  be  charged.  The  resulting 
external  electric  field  will  accelerate  it  to  an  equilibrium  velocity 
such  that  the  Langmuir  condition  is  fulfilled  in  the  frame  of  the 
layer [6,  p.  366]. 

Since  these  criteria  were  derived  using  the  fluid  approximation,  the  pos- 
siblity  that  these  criteria  may  be  relaxed  for  appropriate  velocity  distribu¬ 
tions  remained  to  be  explored.  These  were  the  approaches  of  Kan  and  Lee 
[25]  and  Perkins  and  Sun  [35],  to  be  discussed  below. 


3.5.3  Kan  and  Lee 

Kan  and  Lee  considered  the  same  situation  as  Block  except  that,  instead  of 
including  the  magnetic  field  explicitly,  they  allowed  for  trapped  electrons 


between  the  high  field  side  of  the  DL  and  the  increasing  magnetic  field  of 
the  earth.  For  passing  electrons  the  waterbag  distribution, 

_  j  const  Vi  <  v  <  Vu 
3el  |  0  otherwise, 

was  used,  where  Vuj  —  [(uei  ±  vthei)2  +  but  V}  =  (~)^  for  vei  <  vthe i~ 
similarly  for  passing  ions.  Reflected  ions  and  electrons  were  modeled  with 
Maxwellians.  The  resulting  Bohm  criteria  are 

mv2el  >  2 Tel  +  Tix  Net  <  Nc 
but  veX  >  0  Net  >  Nc 

where  Nc  =  Similarly  A/Vg  >  2Ti0  +  Trt(l  +  ^f).  (Here, 

7  =  2  for  2  degrees  of  freedom.) 

This  result  shows  that  trapped  electrons  can  provide  the  necessary  num¬ 
ber  density  for  electrons  on  the  high  field  side.  However,  we  note  that, 
while  this  approach  indeed  relaxes  the  requirement  on  wei,  Nc  is  a  signifi¬ 
cant  fraction  of  Ne  (since  2e^»  >>  mv2  by  assumption),  and  this  approach 
may  require  an  incident  ion  velocity  significantly  greater  than  the  sound 
speed.  So  while  it  may  be  diminished,  there  still  must  be  a  current. 


3.6  Currentless  Double  Layers 

3.6.1  Perkins  and  Sun 

By  allowing  for  trapped  ions,  Perkins  and  Sun  showed  that  DL  solutions 
exist  which  require  no  current.  In  their  fully  kinetic  treatment  they  chose 
Maxwellian  distributions  for  electrons  and  high  field  ions  and  counter¬ 
streaming  beams  of  trapped  ions  on  the  low  field  side.  The  electron  density 
is  ne  =  n0e~^  where  ip  =  For  equal  ion  and  electron  temperatures, 

r  =  =  1,  the  ion  velocity  distribution  is 

f  _  I  no($r)lje~e  e  >  -A 
3  \  0  c < -A 

2 

This  distribution  requires  the  low-field  particle  kinetic  energy,  =  e  +  tp, 
to  be  greater  than  some  minimum  (t p0  —  A).  A,  then,  determines  how 
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deeply  ions  are  trapped  in  the  potential  well,  ip0.  The  resulting  ion  density 
is 


n0g(ip,A)  =rii  -  /  +  ^)  *de 


[35,  Eqn.  (6)]  and  the  scaled  Poisson’s  equation  (3.7)  becomes 


d2ip 
dX 2 


=  g(ip,  A)  -  e  * 


G&,  A) 


Integrating  by  parts, 

V(<P,A)  =  (||)2  =2j*GU,',AW  >  0 


At  rp0,G(tj>o,  A)  =  0  [35,  Eqn.  8(a)]  for  quasineutrality  (eqn.  3.6), 
while  V(ip o,A)  =  0  [35,  Eqn.  7]  satisfies  Langmuir’s  condition  (eqn.  3.9). 
Eigenvalue  DL  solutions  to  these  equations  exist  for  particular  V’o,  A.  A 
numerical  solution  is  plotted  in  Fig.  3.6,  as  for  Langmuir’s  sheath  of  Figure 
3.4,  the  multivalued  nature  of  <f>{ne  —  rii )  is  to  be  noted.  This  solution  also 
satisfies  Bohm’s  criteria  which  restrict  the  slope  of  the  electron  distribution 
function  on  the  low  density  side  but  are  obviated  by  A  on  the  high  field 
side. 

In  this  presentation  Perkins  and  Sun  make  several  important  observa¬ 
tions.  First,  a  DL  may  exist  in  a  region  where  quasineutrality  imposes  a 
multi-valued  <j> .  Second,  “only”  the  magnitude  not  the  directions,  of  the 
particle  parallel  velocities  are  important.  And  third,  the  particle  distri¬ 
butions  which  support  the  DL  may  be  unstable.  These  observations  have 
important  inference  for  mirrors  and  magnetospheric  potential  theories. 

Their  approach  holds  for  A  a  slowly  varying  (L  A £>)  function  of  x, 
for  example  in  a  mirror  magnetic  field,  and  is  a  generalization  to  analysis 
of  axial  potential  profiles  in  thermal  barrier  cells  [15].  Using  Maxwellian 
distributions  for  electrons  (and  passing  ions),  Cohen  et  al  accounted  for 
trapped  particles  from  energetic  neutral  pumping  by  charge  exchange  with 
a  model  distribution  described  by  a  parameter,  “a”. 


ft  =  n0( 


M 
2t tT 


)?exp[ 


S-ay.  , 

(a  —  1)T 
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♦  =  -e  */  Te 


Figure  3.6:  Electron  and  ion  densities  as  a  function  of  potential  for  r  =  1. 
Curve  a  is  the  electron  density  exp-0.  Curve  b  is  t  3  ion  density 
g(0,A)[Eq.  (6)].  Curve  c  is  the  difference  G(0,A)  and  depicts  regions 
of  positive  and  negative  charge  density.  Dashed  curve  d  would  be  the  ion 
density  if  A  =  0.  It  is  evident  that  the  required  region  of  positive  charge 
density  cannot  exist  for  A  =  0[35,  Fig.  2,  p.  116.] 


Figure  3.7:  Potential  profiles  0(R)  obtained  from  quasi-neutrality  for  vari¬ 
ous  values  of  the  trapped  -particle  filling  parameter  a.  The  portions  of  the 
curves  with  R<  1,  shown  dotted,  apply  when  B>B£^r0oi  as  in  the  dashed 
profile  of  Fig.  1.  The  values  of  a  and  the  corresponding  ratios  p  =  nt/ne 
and  =  —  t<{>/T  evaluated  at  R=3  are:  (A)  a=1.08,  p=0.085,  $m=2.23; 
(B)  a=1.12,  p=0.123,  $m=2.19;  (C)  a=1.15,  p=0.15,  $m=2.15;  (D)  a=1.20, 
p=0.19,  $m=2.08  [15,  Fig.  5,  p.212]. 
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Figure  3.8:  Schematic  view  of  the  experimental  setup.  A  plane  ion  beam 
of  energy  eVb  is  injected  into  the  target  plasma  against  a  biased  magnet 
( Vm  >  Vf,  >  kTe/e).  A  uniform  axial  magnetic  field  is  applied  opposing  the 
magnetic  dipole[44,  p.709]  . 

For  a=0  this  is  a  Maxwellian  while  for  a=l  we  have  a  beam  of  ions.  Plots 
of  quasineutral  solutions,  An($m,  R)  =  0,  are  at  Fig.  3.7.  A  double  layer 
solution  for  appropriate  values  of  the  electric  field  at  its  boundary  is  possible 
where  Langmuir’s  condition  is  satisfied. 

3.6.2  Experiments 

Stenzel,  Ooyama  and  Nakamura  [44]  performed  a  series  of  experiments  in 
which  an  ion  beam  from  a  double  plasma  device  was  injected  toward  a 
dipole  magnetic  field.  As  pictured  in  the  accompanying  sketch  (Fig.  3.8), 
the  beam  was  injected  with  kinetic  energy  (Vj)  by  biasing  one  half  (the 
high  potential  side)  of  a  double  plasma  device  with  respect  to  the  other. 
A  permanent  dipole  magnet  (B  ~  500G)  was  situated  in  the  opposite  half 
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of  the  device  and  a  smaller  axial  magnetic  field  (B  ~  20G)  was  applied 
counter  to  that  of  the  dipole  so  that  the  configuration  consisted  of  open 
and  closed  field  lines  (as  for  the  earth’s  magnetosphere.)  With  magnetic 
potential  Vm  >  Vj  >  e(f>/kT  a  strong  double  layer  was  formed  with  e(f>  ta 
.8 Vb,L  =  <f>/E  ~  lOAp  and  Sn/n  ~  10%. 

In  this  experiment  trapped  particles  were  observed  and  formed  a  large 
percentage  of  the  total  electron  population  ( nt/nj  ~  .48)  [44,  p.715].  The 
role  of  the  magnetic  potential  Vm  was  both  to  totally  reflect  the  ions  and  to 
impart  a  drift  ve  >  vtht  to  the  electrons  trapped  between  the  double  layer 
and  the  magnet.  Of  course,  these  electrons  then  were  susceptible  to  the 
electron  two-stream  (Buneman)  instability.  However,  oscillations  were  not 
recorded  near  the  plasma  frequency,  indicating  a  lack  of  the  ion-acoustic 
instability.  The  authors  attributed  this  to  wave- turbulence  and  an  effective 
T,  >  Te  [44,  p.717]. 

The  experiment  was  performed  with  a  variety  of  plasmas,  the  higher 
mass  ions  being  less  magnetized  .  (The  electrons  are  always  magnetized.) 
Double  layers  were  observed  in  all  cases,  but  the  shape  of  the  DL,  when  the 
ions  were  magnetized,  assumed  the  so-called  inverted  V.  Without  beams 
weak  DLs  were  observed. 

The  double  layer  was  stable.  The  authors  attributed  the  stability  to  a 
momentum  balance  between  the  electrons  accelerated  in  the  double  layer, 
the  reflected  ions,  and  (possibly)  ExB  drifts  in  the  V-shaped  structure. 
This  momentum  balance  is  equivalent  to  Langmuir’s  condition.  Charge 
neutrality  also  plays  a  role  in  that  on  either  side  of  the  double  layer  it  must 
hold [44,  p.717]. 


3.7  Overview 

In  addition  to  learning  useful  methods  for  obtaining  the  local  self-consistent 
solutions  for  potential/particle  systems,  this  chapter  introduced  several  al¬ 
ternative  viewpoints.  On  the  one  hand,  Block’s  (and  Langmuir’s)  work 
demonstrated  the  need  for  satisfying  criteria  on  the  incident  particle  ve¬ 
locity  (Bohm’s)  and  current  density  (Langmuir’s)  for  existence  of  double 
layers.  On  the  other,  theory  by  Perkins  and  Sun,  Cohen,  and  experiments 
argue  for  relaxation  or  removal  of  these  criteria.  In  our  review  we  pointed 
out  where  these  criteria  were  embedded  in  each  of  the  theories  and  exper- 
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AXIAL  POSITION  z  (cm) 


Figure  3.9:  Axial  plasma  potential  profiles  <f>(z )  in  front  of  magnet  with 
bias  Vm.  Without  ion  beam  a  sheath  is  formed,  with  reflected  ion  beam  a 
double  layer  appears  in  front  of  sheath.  Small  magnet[44,  p.712]. 
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imental  conclusions,  and  we  suggest  that  these  differing  points  of  view  are 
really  entirely  consistent.  In  the  next  and  following  chapter  we  shall  fur¬ 
ther  discuss  the  necessity,  or  lack  thereof  for  these  criteria  and  the  ultimate 
harmony  between  competing  points  of  view. 
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Chapter  4 
Analysis 


In  this  chapter,  we  extend  Serizawa  and  Sato  theory  [42]  by  obtaining 
local  solutions  to  the  coupled  Vlasov-Poisson  equations  and  explore  the 
possibility  for  double  layer  solutions  by  applying  the  graphical  and  analytic 
techniques  we  learned  in  the  previous  chapter.  We  remind  the  reader  that 
their  model  is  that  of  a  comoving  neutral  two-species  plasma  injected  along 
the  axis  toward  a  dipole  magnetic  mirror.  However,  in  the  next  chapter 
we  wish  to  use  particle  simulations  to  observe  dynamic  behaviour.  To  best 
accomodate  simulation  limitations,  we  seek  a  simple  model  which  captures 
essential  physical  aspects  we  wish  to  examine  while  permitting  access  to 
analytic  and  computational  solution.  We  devote  this  chapter  to  analysing 
alternative  models  and  comparing  their  implications. 


4.1  Cold  Model 

This  simple  model  is  an  asymptotic  limit  to  the  more  complex  models  we 
discuss.  For  relatively  small  parallel  temperatures,  distributions  introduced 
in  Section  3.6  [35],  [15]  are  approximated  by  beams,  in  which  each  of  the 
particles  have  nearly  the  same  parallel  velocity. 

At  points  other  than  injection  we  expect  the  plasma  to  reach  equilibrium 
according  to  the  Vlasov  equation  (3.1).  It  is  well  known  that  a  function  of 
the  constants  of  motion  must  be  a  steady-state  solution.  In  our  case,  for 
low  collision  frequencies  and  “perpendicular”  frequencies  much  less  them  the 
ion  cyclotron  frequency,  we  take  the  particle  energy,  6  =  +  pB  +  q<f), 
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and  the  magnetic  moment,  (i  =  as  constants  of  motion  and  model  the 
equilibrium  velocity  distributions  accordingly. 

For  a  cold  plasma  with  no  spread  in  perpendicular  or  parallel  velocities 
an  appropriate  distribution  for  each  specie  is 

/  a  S(£  -  £0)8(11  -  Ho) 


or,  rewriting  the  functional  dependence  on  energy  in  terms  of  parallel  ve¬ 
locities  [24,  p.30], 

/  =  A6(v  -  u)6(/j.  -  no)/v 

where  v  =  V||,tt  =  [^(£  -  fiB  -  q<f>)]i  =  -  fiAB  -  qA<j>)]5 ,  and  “0” 

specifies  the  injection  site. 

From  this  distribution  we  calculate  a  density, 

roo  too 

n  —  2n  dv  dv±v±f 
Jo  Jo 


2tt  A 
m 


too  too 

I  dv  dfiBS(v 
Jo  Jo 


)8(h  -  Ho)/v 


2t xAB 


mu 

The  normalization  constant,  A  =  ,  is  evaluated  by  noting  that  n  = 

n0,u  =  v o,  and  B  =  B0  at  injection.  Therefore, 

j  =  nu  =  jo(  — ),  (4.1) 

■Oo 

identical  to  Schmidt’s  analysis  of  a  charge- neutral  homogeneous  plasma 
moving  in  a  strong  magnetic  field  with  weak  curvature  [40,  Eqns.  (6-53), 
p.167]. 

Using  Poisson’s  equation  (3.1(a)),  we  can  evaluate  the  potential  drop  as 
a  function  of  position  (in  terms  of  A B.)  A  particular  solution  results  near 
An  =  n,  —  ne  =  0,  where  the  electric  field,  E,  is  a  maximum  or  minimum. 
Assuming  equal  injected  particle  currents,  jo,  and  magnetic  moments,  /x, 
for  each  specie,  we  obtain  the  electric  potential  in  terms  of  the  magnetic 
field, 

U2  =  u2 

-  h&B  -  eA<£)  =  -(\mvl-H&B  +  eA<t>) 

M  z  m  z 

(\  _  22. 'I 

eA  4>  =  ) - ^4/xA  B 

(1  +  i) 
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where  we  use  capital  letters  to  distinguish  ions  and  lower  case  for  electrons. 
The  maximum  potential  for  this  model  occurs  where  U2  =  u2  =  0  at 

hABmax  =  2^vo~ — 2  ^ 

i  (i  — 

qA<f>MAX  =  2^vl - 2  M  (4-2) 

These  relationships  are  implicit  in  Schmidt’s  analysis  [40,  Eqn.(6-56)]. 

We  observe  that  this  solution  preserves  overall  charge  neutrality  (An  ~ 
0,  Eqn.  (3.6)  )  for  any  <f>.  Quasineutrality  specifies  that  the  total  charge 
contained  even  in  a  double  layer  [45]  or  a  sheath,  when  the  wall  is  consid¬ 
ered,  is  negligible.  In  the  global  sense  then,  quasineutrality  always  holds 
for  the  plasma  [15].  However,  also  assuming  quasineutrality  locally  leads 
to  an  expression  for  the  electric  field 


qE=- 


(i  ~  9MB 
(i  +  »r<fe 


as  a  function  of  the  magnetic  field  gradient.  This  relationship  has  been 
noted  by  many  authors  (e.g.,  [52]  and  [34])  and  results  in  a  potential  de¬ 
picted  in  Figure  4.1.  Also  in  this  figure  is  pictured  a  non-quasineutral 
potential  with  larger  gradients  them  B.  When  the  gradient  in  <f>  is  so  large 
that  the  gradient  in  B  may  be  ignored,  we  may  have  a  double  layer. 

To  further  explore  the  extent  to  which  this  model  preserves  local  quasineu¬ 
trality,  we  develop  a  correction  to  the  quasineutral  solution.  Transforming 
variables  to  V  =  Vi  +  A7  =  7  =  =  (^)3  for  a  dipole  field,  T  = 


L,  and  Ksa  = 


son’s  equation  to 


=  T 


,  simplifies  Pois- 


-V|V,  =  7Ev(>  -  TrT1  +  VJA7  (4.3) 

«  7(^  +  ~)  for  “small”  Vi, 

where  £  =  £0  =  (1  +  §)££,  if  u=U,  and  Z  =  frz  the  axial 

distance  from  the  dipole  center.  For  ^-(z0)  =  0  and  Vj(z0)  =  0, 

■tr  tr  (1  /  ^  \3f21i  .  ( Z  —  Zo)  (Z  —  Z0)n 

Vi  =  Vo{l  +  (y/[-z  sm  L  -  cos  ~L  ]> 
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Figure  4.1:  The  QN  solution  is  E  oc  ^jj-(qA<f>  oc  fiAB).  In  those  cases 
when  the  gradient  of  B  relative  to  E  may  be  ignored,  a  non- quasineutral 
potential,  designated  by  qA</>,  results. 

Far  from  the  injection  point,  Zq  >  Z  >1, 

V,  a  Vo  =  -12(|)! 

remains  “small”. 

The  resulting  electric  field  E  oc  |j  except  near  «  1  where  Equation 
(4.3)  more  resembles  the  ion  sheath  equation  of  section  3.4.1.  This  point 
occurs  near  Z  —  1  where  pff/T  ~  1  and  fj.  may  no  longer  be  treated  as 
a  constant  of  motion.  Although  it  results  in  particles’  ability  to  explore 
a  larger  spatial  extent,  a  potential  doesn’t  necessarily  accelerate  particles. 
The  combined  “potential”  of  the  electric  and  magnetic  fields  slows  both 
electrons  and  ions  in  this  case. 

Using  a  graphical  technique  evident  even  in  Langmuir’s  work,  as  in 
Figure  3.4,  we  obtain  the  quasineutral  solution, 

A  n(  V,  B)  =  0 

Defining  Y  =  ^  =  7  —  1,  in  Fig.  4.2  we  plot  charge  density  A n(V,B) 
for  all  values  V,  Y  <  K.  The  graph  may  be  conveniently  divided  into  four 
sections  separated  by  two  diagonals.  At  the  middle  left  is  the  full  plasma 


Figure  4.2:  This  figure  is  a  charge  density  plot  of  A n(V,B)  in  terms  of 
Y  =  =7—1.  The  light  regions  correspond  to  regions  of  positive 

density  while  the  dark  regions  are  negative.  Along  the  diagonal  extending 
from  the  lower  left  to  the  upper  right  is  the  quasineutral  solution.  The 
white  line  from  the  upper  left  to  lower  right  is  an  ion-rich  (anode)  region 
which  occurs  as  the  ions  overshoot  the  electrons.  No  electrons  are  allowed 
in  the  bottom  quadrant.  Similarly  the  upper  (cathode)  quadrant  permits 
only  electrons  accelerated  by  the  increasing  potential.  The  actual  solution 
of  Poisson’s  equation  is  near  (slightly  below)  the  QN  line  until  it  formally 
becomes  a  sheath  at  the  (anode)  diagonal. 


52 


solution.  Along  the  straight  diagonal  from  the  lower  left  to  the  upper  right 
lies  the  quasineutral  solution  given  by  equation  (3.6).  Above  this  line  the 
solution  is  increasingly  positive  as  the  electrons  accelerate  and  the  ions  slow 
in  the  potential,  V.  Along  the  diagonal  from  the  upper  left  to  the  lower 
right,  is  an  anode  solution  where  the  ions  are  so  dense  that  the  electron 
contribution  may  be  ignored.  In  the  quadrant  below  the  two  diagonals 
electrons  axe  precluded,  as  are  the  ions  above.  To  the  right  neither  are 
allowed. 

As  noted  in  the  previous  chapter,  a  DL  solution  requires  adjacent  layers 
of  oppositely  charged  plasma.  Since  at  fixed  magnetic  field  the  densities 
(4.1)  of  the  two  species  diverge  with  increasing  potential,  these  distributions 
are  incompatible  with  DL  solutions.  However,  as  we  shall  see  in  the  next 
section,  the  assumed  initial  distributions  for  a  particular  situation  may  not 
agree  with  the  steady  state  potential  which  results.  For  example,  in  this 
section  conversion  of  the  cold  electron  velocity  distributions  to  Maxwellian 
leads  instead  to  the  sheath  solutions  of  Section  3.4.1. 


4.2  Finite  Perpendicular  Temperature 

4.2.1  Particle  Currents 


With  finite  perpendicular  temperature  we  must  consider  particles  that  are 
able  to  transit  the  model’s  extent  into  the  loss  cone.  We  choose  a  perpen¬ 
dicular  Maxwellian  injection  distribution 

/  =  /(£,  /*)  =  AS(S  ~  So)  exp  - 
while  leaving  the  parallel  distribution  cold  to  obtain 


J  = 


2t tA 

2  nAT  B 


r  dw  r 

Jo  Jo 


dvlcrc  n 

~^S{S  -So) exp — — 


m“ 


i  fKo~v  JT,  Ko-V-K 

'Jo  dKeXP - Y - 


27 xAT  B  .  Ko-  V. 

where  we  have  used  I\  =  ■ ,  V  =  and  Y  =  ^(=  <fljB  ). 
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At  j  =  jo,  AJ 3  =  0,  B  —  Bq,  and  A (f>  =  0,  giving  A  =  and 

,  ,  B  ^  Ko~V„ 

J  forward  —  JO\  )[1  eXP  (  y  )J’ 

essentially  the  result  produced  in  the  previous  chapter,  Eqn.  4.1,  modified 
by  a  correction  term.  In  the  reverse  direction  we  must  allow  a  loss  cone  for 
those  particles  with  0  <  Ymax  <  Kq  —  Ymax  so  that 

.  /  B  .  Kq  V\fAx  \  — AT0  —  V\ . 

3  re  fleeted  =  M^-)leXP( - 77 - )  “  eXP  ( - 77 - )] 

■d o  Ymax  Y 

for  Ao  V\fAX‘  For  Ao  “C  YmAXi  j  =  j  forward  Jreflected  =  0* 

j  =  Jo(|-)[l  -  exp(-Ao~-^-)]g(Jiro  -  W)  (4.4) 
■do  Ymax 

Requiring  the  steady  state  charge  of  the  system  to  remain  constant,  implies 
A  j  =  ji  —  je  =  0  at  z=0  for  equal  injection  currents  at  z=z0.  Assuming 
equal  injection  velocities, 

eA<f>MAX  Tr  e&<t>MAX 

k«  +  — —  =  K° — fi — 

e&<frMAX  =  2^1(1  -  +  ji’) 

This  is  the  maximum  potential  for  this  model  and  agrees  with  the  cold 
result  (4.2)  for  equal  ion  and  electron  temperatures.  Although  corrections 
for  non-monotonic  <j>  may  be  necessary,  they  were  not  applied  in  the  above 
calculation.  Similarly,  in  the  section  that  follows  we  ignore  corrections 
due  to  the  loss  cone.  However,  as  we  shall  see  in  the  next  chapter,  any 
particular  physical  realization  requires  adjudication  between  these,  possibly 
conflicting,  demands  on  the  potential. 

4.2.2  Particle  Densities 

To  obtain  the  local  solution  for  this  model,  we  begin  by  evaluating  the 
particle  density. 

/oo  roo 

n  =  2n  dv  dv±v±A6(£  —  £0)exp 
Jo  Jo 
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m  A  B 


too  too  „,,1  ,  1  o  .  „ 

J  dv  j  d/jB6(-mv2  —  -mvQ  +  fiAB  +  qA<fi)exp  r 

fOO  bq  imt?—  imv2  —  qA<p .  1  1 

/  dt>exp_AB'  t  )  H{-mvQ  —  -mv2  —  qA(j>) 

Jo  2  2 

(4.5) 


2ttA 
m 

2ttA  B  f00 


„  V  B  T?<  ~  V  \  l 

=  2  n0X0—F( — — — )2 

where  we  have  substituted  for  A  from  above,  H  is  the  Heaviside  function, 
and 

y2  _  Kq  —  V 
Y 

with  Xo  =  X{V  =  0).  The  last  term  in  equation  (4.5), 

F(X )  =  exp  —X2  [  dtexpt2, 

Jo 

is  known  as  Dawson’s  Integral  [1,  p.298]  and  has  the  assymptotic  values 


F(x )  =  xexp—  x2Xl£L0 


2^oo  (-^2)n 


(n  4-  §)«! 


for  small  x  [1,  262]  and 


=  |i  +  £  (1)(3)<^tV)l2n  ~  1)]^for  larEe  *  l4’  p-2561 


n= 1 


In  Figure  4.3  we  have  plotted  n,  and  ne  for  particular  values  of  Y. 
These  establish  the  crossing  points  plotted  in  Figure  4.4,  the  quasineutral 
solutions  for  these  charge  densities.  For  small  (Y,V)  the  solution  is  that 
presented  in  the  previous  section  and  A  <f>  oc  A  B.  For  increasing  A B  an 
additional  solution  appears,  but,  because  n;  increases  with  Y  and  F  reaches 
a  maximum  at  .54,  an  upper  bound  for  Y  exists  at  =  .54 

(treating  n,  as  constant  for  these  values  of  V)  or  y  ~  &-the  point  where  the 
perpendicular  energy  equals  the  total  injected  electron  energy. 

A  solution  also  exists  at  eA<f>  ~  K0  for  all  values  of  Y.  These  solutions 
are  obtained  in  the  (four)  limits  of  large  and  small  x  and  X  bounded  by 
the  two  diagonals  and  lying  in  the  four  quadrants  established  in  Figure  4.2. 
(X,  the  ion  value  for  the  argument  of  Dawson’s  integral,  is  a  maximum  at 
the  lower  left  corner  and  decreases  toward  the  upper  right,  x,  its  electron 
counterpart,  is  a  maximum  at  the  upper  left  corner  and  decreases  toward 
the  lower  right.) 
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Figure  4.3:  Plotted  are  n,-  and  ne  vs  V  for  the  values  Y=.5,  1  and  2.  The 
ion  density  is  that  with  the  peak  to  the  right  and  the  electron  density  is 
the  smaller  peaked  function  to  the  left.  The  points  where  the  two  densities 
cross  represent  quasineutral  solutions. 
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Figure  4.4:  Plotted  here  are  the  contours  An  =  0  for  both  electron  and  ion 
densities  proportional  to  F(X). 
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Between  the  two  quasineutral  solutions  n,  >  ne  and  a  sheath  must  ex¬ 
ist.  To  connect  these  solutions  requires  passage  through  a  region  of  strong 
electric  field  gradient.  As  noted  by  Stenzel  and  Ooyama  [44,  p.714],  accel¬ 
eration  of  electrons  in  this  region  may  lead  to  a  strong  electron  two-stream 
instability  and  their  rapid  thermalization.  Designating  this  point  as  “1”, 
we  may  replace  the  Dawsonian  electron  distribution  with  a  Maxwellian 

n  =  exp  ^ 

where  n\  =  2hq^XqF{ Kaf v- ) \ Zl  and  Te  represents  a  spread  induced  by 
thermalization,  expected  to  be  on  the  order  of  the  accelerating  potential. 

Plotted  in  Fig.  4.5  is  the  contour  An  =  0  for  K  =  25,  Te  =  15  and  z\ 
at  Y=l.  Near  z\  the  potential  V  is  double- valued,  allowing  a  DL  solution 
when  Langmuir’s  condition,  eqn.  (3.9),  is  also  satisfied.  We  remark  that 
the  shapes  of  these  model  distributions  lead  to  double  layer  solutions,  as 
Perkins  and  Sun  [35]  postulated  and  we  demonstrate  in  the  next  chapter. 


4.2.3  Langmuir’s  Condition 

Holding  B  constant,  an  exact  integral  of  F(X)  over  <f>  is 

e  (*  F(X)d<j>/T  =  [V  dV  exp  - X 2  [*  dt  exp  t 2 
Jo  Jo  Jo 

=  f  dV  exp  —  V  f  dt  exp  t 2 
Jo  Y  Jo 

=  Y{F(X)-X\\V0, 

This  result  is  useful  for  establishing  Langmuir’s  condition: 

qJ(ni-ne)d4>  =  Ti{2n0^-QXo[F(X)-X}\l}  -  Te(n  -  n0) 

Ti(n  -  n0)  -  Te(n  -  n0)  -  (X  -  X0)(2n0X0~)  =  0 

-oo 

AT  An  -  AX(2n0X0~)  =  0 

Bq 

Together  with  quasineutrality,  this  condition  yields  a  double  layer  solu¬ 
tion,  visually  evident  in  Figure  4.5  when  the  integrated  particle  densities 
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Maxwellian  Electrons  an^Dawsonhm  Ions^^Th  C°n*'°Ur  An  =  0  fo 
ble- valued  in  V,  and  a  DL  oni,  *•  1  The  soIutlon  is  obviously  dou 

satisfied.  DL  Solut,on  -here  Langmuir’s  condition  " 


over  the  potential  enclosed  by  two  areas  on  either  side  of  a  line  drawn  per¬ 
pendicular  to  the  Y  axis  and  the  quasineutral  solution  are  equal.  Since 
by  assumption  AX  <  0  and  An  >  0,  Te  >  T{  =<  fii  >  A B  is  required 
to  satisfy  Langmuir’s  condition,  consistent  with  streaming  instabilities[47]. 
While  this  distribution  will  prove  valuable  in  our  particle  simulations,  in¬ 
stabilities  will  permit  observation  of  the  more  general  model  we  discuss 
below. 


4.3  Generalization  to  Drifting  Maxwellian 

4.3.1  Current  Density 

We  generalize  our  results  to  a  drifting  Maxwellian  by  integrating  over  the 
distribution 

h(u)  =  xexP-!n(;;-u)2 
with  the  boundary  condition 

f<x> 

jo  =  no  /  h(u)udu  =  n0  <  u  > 

Jo 

for  both  species.  Using  Eqn.  (4.4),  we  evaluate  the  particle  current  density 
as 


B  f°° 

j  —  no(ir)  duuh(u)[  1-exp 

&0  Jum,n 


K  -  V, 


B  t°°  K  —  V  rumin 

n0 {~B~)[<U>-  duuh(u)e\ p - — - /  duuh(u)(  1  -  exp - 

B  o  Jo  Y  Jo 


K-V, 


Jo  Jo 


where 


—  | 


(2 ej>/m)3  for  ions 


for  electrons 


and  we  redefine  T  =  T\\  and  /i  =<  n  >.  K,  V  and  Y  have  their  former  values 
multiplied  by  the  quantity  For  example,  V  =  =  Kmin. 

We  recognize  j0  as  the  cold  particle  current  density,  Eqn.  4.1. 
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For  Kq  —  V  1,<  u  >«  ^£.y/K^f£°dAxexp—Ax2,  and  we  may 
approximate  A j 
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where  B  s  =  dfe-  Similarly, 
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The  ratio  of  the  last  two  terms  is  |^4|  =  y  f°r  large  Y. 

Requiring  In  x  =  0  for  charge  conservation,  we  set 

Aj,  A  j  ,  jj 

x  =  —  =  +  X7'  =  xo  +  xi 

Aje  Aje  Aje 

Since  lnx(Vo  +  Vi)  «  lnx0(Vo)  +  +  Eaf^kYiJ  =  0,  to  lowest 

order  lnx0(VJ))  =  0  and 
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In  the  limit  of  large  injection  energy  to  parallel  temperature  ratio,  this 
assymptotic  evaluation  of  the  potential  has  resulted  in  an  analytic  expres¬ 
sion  which  agrees  exactly  with  the  “semi-empirical”  result  obtained  by  Ser- 
izawa  and  Sato  [42,  Eqn.  7]  but  provides  corrections  which  account  for 
biMaxwellian  temperature  distributions. 

4.3.2  Particle  Density 

We  compare  the  density  for  this  distribution  to  that  of  the  previous  sections. 
Integrating  over  (4.5) 

n  =  2n0—  j^_F{ — — — )?  exp  -(z  -  \jK0)2dz(y )?/  jf  exp  -(z  -  yjK0)2dz 

The  integral  can  be  evaluated  by  a  variation  on  the  method  of  steepest 
descent[4].  We  begin  by  changing  variables  to  t  =  z/\/V,r  =  Noting 

that  F(X)  may  be  written 

F(X)  —  exp  —  X2  f  dxex px2  =  X  /  dyexp —(l  —  y2)X2, 

Jo  Jo 

n(t,y)  =  2n0^--jL=:J'  dt  jf  dy[y(t2~l)}hexp -[(1  -  y2)y(t2  -  1)  +  V(t  -  r)2] 

Setting  4>(t,y)  =  (1  —  y2)(t2  —  1)  +  Y(t  —  r)2,  we  note 

V0t,v  =  0  at  {<0,yo}  =  {1, [1  +  F(1  -  r)jj}  and  {t0,3/o}  =  {j^O},  where 
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These  two  points  correspond  to  the  combined  contribution  from  the  two 
peaks  in  F  and  the  Maxwellian. 


For  t0  <  1  n  = 
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Letting  A y  =  mAt,  we  integrate  along  a  path  of  steepest  descent  in  the 
Ay,  At  plane  for  m  =  —  ^  yielding 
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with  I(x)  =  i[(l  +  x)2  -  l]t,  xmax  =  — m[(l  -  y0)  +  P],  and 


X 


min  — 


( 


0  -  (1  -  yo)  < 
—m(—y0  +  P)yo  < 


P 

P 


<  yo 

<  oo 


so  that  n  decays  exponentially  from  r=l. 
Similarly  for  t0  >  1, 


j|p)!  exp -V(t0  -  r fF[~(t\  -  l)]i, 


which  corresponds  to  the  original  distribution  shifted  in  Y  and  diminished 
by  an  exponential  factor.  Finally,  for  to  ^  1  and  y0  ~  0  the  contribution 
from  the  two  peaks  coalesce  giving 


n  «  2n0 
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A  slowly  converging  function  yielding  a  broad  peak. 

These  three  equations  correspond  with  the  original  quadrants  in  Figure 
4.2.  The  last  essentially  divides  the  plane  along  the  quasineutral  line,  above 
and  below  are  regions  of  small  and  large  t0  respectively,  apparent  as  the 
high  V  and  low  V  areas  in  Figure  4.6,  where  we  have  plotted  the  densities 
for  Y=1  and  4.  We  see  that,  as  expected,  the  peak  is  spread.  Although  the 
inflection  point  is  lowered,  it  permits  the  ion  density  to  cross  a  constantly 
increasing  ne  at  three  points  in  V.  The  shift  in  the  peak  lowers  the  range 
of  n  values,  but  not  necessarily  the  magnitude  of  double  layer  potentials. 
Just  as  in  the  previous  section,  opportunity  for  a  DL  solution  diminishes 
with  increasing  B.  We  shall  observe  the  consequences  of  this  analysis  in  the 
next  chapter. 


63 


Figure  4.6:  In  this  figure  are  plotted  ion  distribution  functions  where  in¬ 
jected  velocities  are  drawn  from  a  Drifting  Maxwellian.  Comparison  to  the 
F  distribution  of  the  previous  section  shows  that  the  peak  in  this  distribu¬ 
tion  is  less  pronounced.  Similarly  the  inflection  point  has  shifted  to  the  left 
(toward  injection).  However,  presence  of  the  inflection  point  together  with 
a  sufficiently  large  electron  temperatures  should  allow  DL  solutions  in  the 
cases  we  consider. 
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4.4  Summary 

The  analysis  of  the  potentials  we  developed  in  these  sections  contrasts  be¬ 
tween  the  global  requirement  which  maintains  overall  charge  neutrality  and 
the  local  one  which  permits  double  layer  solutions.  The  latter  potential  is 
dependent  on  the  self  consistent  constituent  particle  distributions  both  in 
slope  (in  terms  of  Bohm  criteria)  and  in  integrated  charge  density  over  the 
potential  (in  terms  of  Langmuir’s  condition.) 

The  results  of  this  section  showed  that,  in  limited  circumstances  of  the 
cold  one-dimensional  plasma  in  presence  of  a  magnetic  field,  the  only  po¬ 
tential  solution  is  one  in  which  the  plasma  must  either  reach  quasineutrality 
or  leave  a  region  unexplored.  First  order  corrections  show  explicitly  that 
large  scale  potential  drops  do  not  result  in  double  layers  or  even  particle 
acceleration.  Thus,  the  global  and  the  local  requirements  are  equivalent. 

With  the  added  freedom  to  redistribute  energy  between  parallel  and 
perpendicular  velocities  the  plasma  may  explore  beyond  boundaries  for  the 
cold  plasma.  Indeed,  the  plasma  has  access  to  an  infinitude  of  quasineutral 
solutions.  The  alternative  it  chooses  is  that  required  for  pressure  balance 
as  demanded  by  Langmuir’s  condition.  However,  the  conditions  for  double 
layers  are  stringent,  and  for  large  mass  ratios  a  double  layer  solution  is 
difficult  to  attain  with  Dawsonian  electrons. 

When  instabilities  alter  particle  distributions,  as  we  allowed  in  the  case 
of  cold  parallel  electrons,  they  may  present  opportunities  for  double  layer 
solutions.  Changing  the  electron  distribution  in  Section  4.2  to  that  of  a 
Maxwellian  immediately  leads  to  double  layer  solutions  pictured  in  Figure 
4.5.  These  distributions  and  results  agree  with  those  first  proposed  by 
Perkins  and  Sun  [35].  Addition  of  parallel  temperatures  softens  distribution 
shapes,  but  still  permits  double  layer  solutions  to  exist.  Indeed,  apparent 
from  graphical  methods,  local  perturbations  may  even  provide  opportunity 
for  multiple  closely  spaced  double  layers. 

In  each  case  we  demonstrated  Langmuir’s  condition  in  the  role  of  pres¬ 
sure  balance  and  couched  Bohm’s  criteria  as  a  requirement  on  the  slope 
of  the  distributions  functions.  This  view  provides  different  emphasis  from 
but  demonstrates  the  consistency  between  frameworks  established  by  other 
authors. 

Our  global  analysis  confirms  Serizawa  and  Sato  [42].  We  demonstrate 
that  their  solution  varies  little  between  cold  (Eqn.  4.2)  and  warm  (Eqn. 
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4.6)  paxticles  except  for  a  multiplicative  factor  of  order  unity,  attesting 
to  a  weak  dependence  on  the  mirror  ratio  Tj|  brings.  Additionally,  the 
physical  distribution  of  Serizawa  and  Sato  in  section  4.3  should  permit 
double  layer  solutions.  As  demonstrated,  the  model  distribution  of  section 
4.2  well  approximates  both  the  global  and  local  results  expected  with  their 
distributions.  However,  to  observe  the  time- dependent  behaviour  of  these 
scenarios,  we  resort  to  particle  simulation  as  presented  in  the  next  chapter. 

Actual  solutions  depend  greatly  on  not  their  initial  values  but  on  the 
time-dependent  behaviour  of  the  resulting  self-consistent  particle  distribu¬ 
tions  actually  explored.  To  investigate  these,  it  is  necessary  to  resort  to 
particle  simulation  techniques.  In  Chapter  5  we  investigate  a  simulation 
equivalent  to  the  cold  7j|,  warm  T±  distribution  of  section  4.2.  As  pointed 
out  above,  the  adequacy  of  this  model  was  demonstrated  by  its  accessibility 
to  double  layer  solutions  and  the  close  resemblance  it  bears  to  the  expected 
analytical  results  when  a  thermal  parallel  distribution  is  added.  In  the  sub¬ 
sequent  chapter  we  extend  our  one-dimensional  results  to  preliminary  ones 
in  two  dimensions. 
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Chapter  5 


One  Dimensional  Particle 
Simulation 


In  the  previous  chapter  we  explored  the  static  local  and  global  solutions 
to  Poisson- Vlasov  equations.  This,  as  well  as  previous  theory,  lacks  from 
the  ability  to  observe  time  dependent  behaviour-instabilities,  stationarity, 
and  accessibility.  While  experiments  allow  observation  of  these  phenomena, 
often  they  suffer  from  constraints  of  physical  accessibility,  range  of  param¬ 
eter  exploration  and  cost.  In  this  chapter  we  compare  the  solutions  in  the 
previous  chapter  to  those  using  particle  simulation. 

Previous  simulations  were  limited  in  variety  of  particle  and  field  bound¬ 
ary  conditions  allowed.  Our  model  advances  beyond  these  to  allow  injection 
with  floating  potentials  in  a  bounded  simulation  so  that  the  double  layer 
isn’t  predisposed  by  initial  conditions.  While  the  electric  field  is  calcu¬ 
lated  in  one  dimension  only,  thus  ID,  the  magnetic  field  is  a  realistic  one 
required  for  observing  effects  of  the  magnetic  moment.  Once  injected,  par¬ 
ticles  are  free  to  explore  all  of  position  and  velocity  space  allowed  by  the 
one-dimensional  potential  and  the  three-dimensional  magnetic  field. 


5.1  Physical  Model 

The  physical  model  is  that  of  auroral  field  lines  with  one  end  tied  near 
one  of  the  earth’s  magnetic  poles  and  the  opposite  embedded  near  the 
reconnection  zone  far  from  the  earth,  but  this  model  may  apply  equally  to 
other  “mirror”  configurations. 
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Figure  5.1:  A  projection  of  the  simulation  region.  Particles  are  injected  into 
a  box  at  the  right  boundary  (z=L)  where  symmetry  requires  the  electric 
field  to  be  zero.  The  axis  of  the  simulation  is  that  of  the  dipole  magnetic 
field.  The  potent’ J  on  the  left  boundary  (z=0)  is  defined  as  <f>= 0. 

5.1.1  Particles 

Currentless  plasma  enters  the  modeled  region  from  the  reconnection  zone 
streaming  towards  the  earth  with  injection  velocities  vj  ~  lOOOfcm/s,  yield¬ 
ing  ion  kinetic  energies  K  ~  \0keV.  Temperatures  are  Te  <  T,  ~  1  keV. 
Both  the  density  and  magnetic  field  are  functions  of  distance  from  the 
earth,  R.  Extrapolating  va  =  500 km /a  [38,  p.  7180]  in  the  tail  region  of 
the  magnetosphere  where  R  ~  10 Re  to  the  ionosphere  where  R  ~  Re  and 
using  va  oc  and  n  oc  B  yield  densities  n  ~  100cm-3. 

In  this  region  particle  and  field  scale  lengths  are 

~  Pe  pi  R,  \d  <  l  £  R, 

where  R  is  the  B  field  scale  length  and  l  is  that  expected  for  E.  Similarly, 
the  frequency  ordering  is 

<  Upi  ~  U)  <C  ~  ^pe> 

fl  the  gyrofrequency.  Particle  collisions  may  be  ignored. 


5.1.2  Electric  Field 

For  these  parameters,  {3  <C  1.  In  this  regime  electrostatic  (ES)  and  electro¬ 
magnetic  (EM)  waves  are  decoupled  [46,  p.  224]  and  [28,  p.  427],  and  we 
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may  choose  an  ES  solution  to  Maxwell’s  equations.  The  boundary  condi¬ 
tions  for  the  electric  potential  are  chosen  as  (f)  =  0  at  z=0  (the  ionosphere  or 
mirror  throat)  and  E=0  at  z=L  (the  magnetosphere  or  mirror  center.)  As¬ 
suming  that  the  same  physical  processes  drive  phenomena  at  both  the  North 
and  South  poles  requires  these  potentials  and  charge  distributions  to  mir¬ 
ror  one  another  about  the  equatorial  plane  of  the  earth.  While  disallowing 
odd  solutions  to  Poisson’s  equation,  this  assumption  leads  to  the  condition 
E(L)=0.  The  condition  at  z=0  derives  from  the  presence  of  high  density 
regions  of  ionized  gas  near  the  earth’s  poles.  Taking  the  poles  as  electric 
ground  while  assuming  local  charge  neutrality  requires  E=0  throughout  the 
regions  and  </i>(0)  =  0.  (For  discussion  of  these  boundary  conditions  see  [15] 
and  [48].) 


5.1.3  Magnetic  Field 

We  use  the  exact  dipole  field  of  Section  2.2.1.  However,  taking  the  electric 
field  scale  as  K/l  with  /  <C  R  for  a  double  layer  and  mv 2  ~  T, 


mv2/R 

qE 


«  1, 


permitting  us  to  ignore  curvature  drifts  by  centering  the  model  and  aligning 
the  direction  of  the  electric  field  along  the  magnetic  field  axis.  The  value 
of  the  field  is  chosen  so  that  at  the  model’s  center  the  magnitude  of  the 
field  is  equal  to  that  at  R  =  2.5Re  directly  above  the  poles.  However,  to 
reduce  the  size  of  the  simulation  region,  the  magnetic  field  scale  length  is 
reduced  while  maintaining  the  above  relative  parameters.  Bounding  the 
computational  model  are  injection  and  loss  planes  perpendicular  to  the 
magnetic  field  axis. 


5.2  The  Computational  Model 

We  primarily  choose  parameter  values  for  our  simulation  to  maintain  fre¬ 
quency  and  length  scalings.  In  part,  these  values  mirror  those  of  the  physi¬ 
cal  model  (Figure  5.1.)  However,  a  particle  simulation  is  limited  by  numbers 
of  particles  and  simulation  size  and  duration.  Therefore,  use  of  this  tool 
requires  careful  abstraction  of  both  physical  processes  and  parameters. 


5.2.1  Electric  Field 

Electric  Field  Grid 

To  calculate  the  electric  field  we  divide  the  computational  model  at  equi- 
spaced  grid  points  defined  by  intersection  of  intermediate  planes  perpen¬ 
dicular  to  the  magnetic  axis.  At  each  time  step  particles  are  assigned  to 
grid  points  centered  between  planes  to  obtain  charge  density.  The  electric 
field  is  solved  and  the  particles  are  moved  according  to  the  electric  force 
interpolated  to  their  positions. 

The  grid  length,  A,  is  chosen  as  Ad.  ID  simulation  theory  for  no  (or 
constant)  magnetic  field  has  demonstrated  the  effectiveness  of  this  choice 
[23].  In  the  presence  of  magnetic  field  gradients,  however,  the  density  is 
not  uniform  and  one  may  propose  the  smallest  “Ad”  to  capture  fine  scale 
behaviour  at  the  highest  densities.  We  have  chosen  A d  at  the  injection 
point. 


Charge  Sharing/Force  Interpolation 

In  our  ID  simulation  we  use  an  interpolation  scheme  which  assigns  particle 
charges  to  grid  points,  xg, 
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and  interpolates  particle  force  from  electric  fields  at  gridpoints,  Eg, 


Fi  —  TlgSigEg 

Because  electric  fields  are  calculated  in  the  axial  direction  only,  the  same 
equations  are  obtained  by  considering  particles  to  be  constant  density  slabs 
of  infinite  cross  section  but  finite  thickness,  A,  rnd  basing  charge  shar¬ 
ing/force  interpolation  on  the  particle  volume  in  each  cell.  This  scheme  is 
known  as  cloud-in-cell  (CIC)  weighting  [5,  p.  311]. 


Poisson  Solver 

The  charge  sharing  scheme  above  results  in  a  grid  charge  density 

Pg  — 
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where  the  summation  is  over  individual  particles,  i. 

Direct  integration  of  Poisson’s  equation  3.1(a)  in  ID, 

Eg  —  Eg+l  =  i  A , 

is  fast  and  easily  implemented.  However,  especially  for  diagnostics,  ad¬ 
vantages  accrue  by  using  a  Fourier  solver.  Some  authors  have  advocated 
straightforward  application  of  continuous  Fourier  transforms  to  the  dis¬ 
crete  case  [17].  In  general,  these  methods  don’t  account  for  effects  of  alias¬ 
ing,  which  may  predominate  in  bounded  simulations.  After  correcting  for 
these  aliases,  using  techniques  developed  in  Appendix  B,  their  effects  are 
ameliorated  and  the  Fourier  transform  results  closely  agree  with  those  ob¬ 
tained  by  direct  integration.  Thus,  direct  integration  was  used  for  the  ID 
simulation  while  post-processing  diagnostics  incorporated  corrected  Fourier 
transforms. 


5.2.2  Particles 

Particle  Pusher 


While  taking  advantage  of  axial  symmetry,  particle  acceleration  is  calcu¬ 
lated  with  a  full  pusher  or,  for  electrons,  a  guiding  center  (gc)  pusher.  In 
this  chapter  we  primarily  report  results  from  the  full  pusher, 
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This  equation,  along  with  its  position  counterpart, 


xi  —  Xo 
At 


constitutes  the  “Leap-Frog”  scheme  (Eqn.  A.l). 

In  the  next  chapter  we  report  results  in  2D  using  a  gc  pusher  (Eqn.  A. 3). 
The  full  and  gc  pushers  yield  similar  results  in  ID.  In  practice,  however, 
the  charge  to  mass  ratio  (q/m)  is  chosen  to  retain  the  appropriate  scale 
lengths  while  relaxing  the  required  number  of  time  steps  to  observe  the 
desired  physical  phenomena.  At  injection  care  is  taken  to  retard  particle 
positions  by  a  half  time  step  in  adherence  to  the  leap-frog  scheme. 
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Particle  Velocity  Distributions 

Particles  are  injected  at  the  low-field  side  of  our  cylindrical  box  of  length 
L  =  N\p.  To  simulate  a  warm  plasma  accelerated  by  external  processes, 
the  injected  particles  would  be  drawn  ideally  from  a  drifting  Maxwellian 
velocity  distribution.  However,  previous  authors  demonstrated  “energy  in¬ 
stabilities”  attributed  to  current  flow  through  a  boundary  at  floating  poten¬ 
tial  [48].  Because  we  expected  a  time-dependent  potential  at  the  injection 
boundary,  we  sought  a  computational  model  which  retained  features  of  par¬ 
ticle  mirroring  and  drift  while  allowing  simple  and  current-free  injection. 
For  both  ions  and  electrons  we  use  a  cold  (delta-function)  distribution, 

fo  =  An06(v  -  a)  exp 

the  injection  probability  is  simply  P=1  for  v  =  U||  =  u  and  0  otherwise. 
This  approach  permits  a  simple  numerical  algorithm  for  particle  injection 
while  avoiding  these  “energy  instabilities.”  As  presented  in  Section  4.2,  this 
distribution  allows  straightforward  comparison  to  results  expected  from  a 
more  general  drifting  Maxwellian.  Indeed,  while  this  distribution  results  in 
simulations  consistent  with  our  analyses,  thermal  processes  permit  direct 
extension  of  our  results  to  more  general  situations  than  this  simple  approach 
might  suggest.  These  aspects  will  be  discussed  below. 

Particle  Injection 

The  number  of  particles  penetrating  a  boundary  of  cross  section,  A,  during 
each  time  step,  At,  is 

J  dA  ■  J  fvdvAt  =  TAfA 

The  flux  for  an  azimuthally  symmetric  distribution  along  the  magnetic  field 
axis, 

r  =  7 T  J  dvlj  <ty|/(U||,Vj>||, 

is  weighted  by  its  parallel  velocity. 

At  any  particular  time  step  the  probability  for  injection  with  perpen¬ 
dicular  velocity  of  magnitude  greater  than  =  |tTx  |  is 

roo  too  too 

P(v  1)=  .  dv\ g(vi_)/  /  dv±g(v±)  where  g  =  /  dv\\f(vn,v±) 

Jv,  JO  J—oo 
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Noting  that  P(0)=1,  to  distribute  velocities  for  a  perpendicular  Maxwellian 
distribution  (g  oc  exp  —mv±/2T),  we  simply  choose 

v±  =  ~~  In  P(x) 

m 

where  P(x)  is  a  uniform  probability  distribution  and  x  is  a  pseudo-random 
number  distributed  on  the  interval  x  e  [0,1]. 

In  our  simulation  the  particle  injection  rate  is  N /  At  =  n0voA.  This 
quantity  is  fixed  by  choice  of  the  three  parameters,  uo,vq,  and  At.  To 
inject  discrete  particles,  one  may  either  choose  these  parameters  for  integer 
injection  number,  N,  or  he  must  account  for  residual  particles  left  behind 
at  each  time  step.  To  maintain  a  currentless  boundary,  the  actual  injection 
was  most  often  performed  by  placing  electron  guiding  centers  on  top  of  ions’ 
and  injecting  a  fixed  number  of  these  pairs  with  identical  parallel  velocities 
at  each  time  step. 


Particle  Boundary  Conditions 


Particles  that  exit  at  z=L  are  mirrored  to  the  position  2(L-z)  with  perpen¬ 
dicular  velocities  identical  to,  but  parallel  velocities  opposite,  the  previous 
time  step.  Various  schemes  were  used  at  z=0. 

Replacing  particles  at  z=0  without  considering  their  attributes,  such  as 
“recycling”  particles  by  exchanging  those  lost  at  z=0  with  replacements  at 
z=L,  generally  do  not  conserve  energy.  Therefore,  particles  lost  at  z=0  were 
most  often  simply  allowed  to  escape,  requiring  the  potential  to  maintain 
overall  system  charge  neutrality.  That  is,  a  sheath  was  formed. 

One  alternative  was  to  replace  hot  plasma,  lost  during  time  scales  of 
interest,  with  backscattered  ionospheric  particles  cooled  by  collisions  with 
neutrals.  Steady  state  requires  that  the  total  flux  of  hot  charge  balances 
that  of  cold.  Approximating  v  ss 


r  cold 


l^cold 


-r  hot 


(nv)hot 

VCold 


This  technique  was  used  to  study  the  effects  of  these  fluxes  on  simulation 
results  and  agrees  with  results  we  report  below  except,  as  noted,  the  ap¬ 
proximation  Ar  =  0  ameliorates  the  need  for  a  sheath. 
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5.3  Simulation  Results 

5.3.1  Simulation  Parameters 

We  present  results  from  a  simulation  with  Ng  =  L/Xd  —  1024,  M/m  =  25, 
uJpeAt  =  0.25,  and  mirror  ratio  7  ~  25  with  B  at  the  center  point  of  the 
model  equal  to  that  at  R  —  2.5 Re  and  an  injection  density  no  =  lOOcm-3 
Perpendicular  velocities  are  normally  distributed  with  Tt  =  2 Te  =  1  keV. 
The  ions  and  electrons  share  a  common  injection  velocity  V{nj  =  .6y/2vthe- 
These  values  yield  upe  ~  Q,e  ~  105,  A^  ~  pe  ~  10m,  and  p,  ~  100m.  Even 
in  the  absence  of  particle  collisions  the  ion  bounce  time  and  loss  cone  size 
determine  the  extent  to  which  the  ion  velocity  distribution  is  characterized 
by  the  injected  beam.  A  crude  estimate  of  the  bounce  time  is  rj,  ~  1000/u;pe 
which  allows  sufficient  time  for  the  bulk  of  the  injected  ions  to  be  reflected 
by  the  magnetic  field,  but  only  few  ions  may  complete  an  “orbit”  in  z. 

5.3.2  Chronology 

In  Figure  5.2  are  plots  of  velocity  space  vs  axial  position  at  times  upet  = 
2100  and  3300  for  a  charge  neutral  plasma  injected  at  z=L.  The  instanta¬ 
neous  relative  charge  density  and  electric  field  are  also  plotted  and  can  be 
compared  to  the  simulation’s  progress  in  velocity  space. 

As  the  plasma  drifts  into  the  dipole  magnetic  field,  a  strong  electric 
field  is  evident  by  the  acceleration  of  electrons  and  slowing  of  ions  at  z  w 
700 A/)  for  ujpet  =  2100  and  z  ~  750 A#  for  upet  =  3300.  Acceleration  is 
apparent  in  the  growth  of  phase  space  and  increase  in  average  velocity  to 
<  ve  >«  1.25— 1.5vt/,e  at  u)pet  =  2100  and  <  ve  >«  2-Ou^e  at  upet  =  3300  for 
electrons  and  by  their  contraction  and  decrease  for  ions.  The  instantaneous 
electric  field  and  charge  densities  are  maximum  near  these  points,  indicating 
a  potential  drop. 

Electrons  are  rapidly  scattered  in  parallel  velocity  while  ions  slowly  fill 
velocity  space  behind  the  potential.  The  width  of  the  electron  distributions 
implies  that  their  effective  parallel  temperature  is  much  greater  than  the 
injection  temperature.  Underlying  the  velocity  space  plots  are  oscillations 
Ave  ~  vthe  with  scale  lengths  A  ~  10-sAp  for  both  ions  and  electrons.  One 
also  notes  that  the  potential  moves  toward  the  “equator”  at  a  velocity 
vdl  ~  (100Ao)/(1200/u;p«)  ~  .lvi&e  <  vthi  and  that  its  movement  coincides 
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Figure  5.2:  Plotted  in  these  two  sets  of  panels  are  phase  space  ( vz  —  z) 
(top)  and  average  velocities  (bottom)  for  electrons  (left)  and  ions  (center)  in 
units  of  electron  thermal  velocity,  plots  of  the  electric  field  acceleration  per 
time  step  (change  in  parallel  velocity)  normalized  to  the  thermal  velocity 
(top  right),  and  instantaneous  charge  density  normalized  to  the  injected 
particle  density  averaged  over  several  cells  (bottom  right).  The  spatial  grid 
is  graduated  in  units  of  s  =  z/Xq  =  100.  For  this  simulation  particles  are 
injected  at  s=1024  on  the  right  with  the  “ionosphere”  on  the  left  at  s=0. 
The  plots  are  for  two  different  times  upet  =  2100  (above)  and  3300  (below). 


with  the  filling  of  ion  velocity  space. 


5.3.3  Velocity  Distributions 

Shown  in  Figure  5.3  are  electron  and  ion  velocity  distributions  at  times 
u>pet  =  2100  (top)  and  3300  (bottom)  at  eight  equally  spaced  grid  points 
from  s  =  z/\d  =  64  (upper  left)  to  960  (lower  right).  For  large  z  two  ion 
beams  are  clearly  present.  The  peaks  of  these  beams  converge  until,  beyond 
the  potential  drop  of  Figure  5.2,  they  are  indistinguishable.  The  assymetry 
in  the  electron  distribution  at  large  z  indicates  velocity  components  in  both, 
but  particularly  the  forward  (injection),  directions.  For  intermediate  z  val¬ 
ues  the  effect  of  the  electric  field  is  evident  where  the  electron  distribution 
spreads  and  flattens.  At  smaller  z  the  distribution  spreads  to  greater  and 
greater  widths  but  then  again  narrows. 

Accounting  for  differences  in  mass,  the  spread  of  the  two  species’  veloc¬ 
ities  are  comparable  at  large  and  very  small  z,  indicating  roughly  the  same 
temperatures.  However,  near  the  potential  drop  the  ions  remain  beams  and 
their  velocities  spread  relatively  less  than  those  of  the  thermal  electrons. 
These  observations  are  borne  out  in  the  next  two  figures  (5.4  and  5.5) 
obtained  for  a  gc  simulation  with  parameters  identical  to  this  simulation 
except  for  jh  —  1  at  injection. 

5.3.4  Stringer  Plots 

Figures  5.4  and  5.5,  as  those  presented  by  Stringer [47],  display  velocity 
differences  and  temperature  ratios  as  functions  of  position  for  each  pair 
of  constituents  in  counterstreaming  plasmas.  One  set,  Figure  5.4,  por¬ 
trays  differential  (current-carrying)  flows  between  ions  and  electrons.  The 
other  set,  Figure  5.5,  displays  differential  ion-ion,  electron-electron  and  ion- 
electron  flows  between  counterstreaming  neutral  beams.  At  right  of  each 
set  of  differential  velocity  and  temperature  ratio  plots  are  their  point  values 
for  each  of  several  positions  along  the  axis.  Comparison  of  these  plots  to 
regions  of  instability,  while  considering  differences  in  ion  mass  and  a  factor 
of  two  difference  in  like  species  counterstreaming  velocities,  were  intended 
to  identify  possible  regions  of  instability. 

Stringer’s  plots  were  based  on  drifting  Maxwellian  distributions  and  a 
realistic  mass  ratio  (1836).  But,  as  discussed  in  Section  5.4.3  below,  our 
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Figure  5.3:  Relative  densities  of  electrons  and  ions  in  velocity  space  grad¬ 
uated  in  units  of  vj he-  77 
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Figure  5.4:  These  are  current  carrying  flows  at  a >pet  —  2600.  The  top  figure 
displays  positive  flowing  ions  and  negative  flowing  electrons.  The  bottom 
is  the  opposite  comparison.  Differential  velocities  are  in  units  of  vthe  while 
the  position  is  graduated  in  units  of  100s  of  Debye  lengths. 
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Figure  5.5:  The  same  comparisons  as  the  previous  figure’s  but  for  counter¬ 
streaming  neutral  beams.  The  top  left  panel  represents  counterstreaming 
ions  (positive  velocities)  and  electrons  (negative  velocities).  Because  the  e-e 
wave  is  unstable  at  any  temperature  for  Av  >  2.6vthe,  it  is  omitted  from 
the  adjoining  two  plots.  The  two  bottom  rows  represent  the  negative  (top) 
and  positive  (bottom)  ion  beams  interacting  with  the  combined  electron 
beams. 
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simulation  certainly  leads  to  other  nearby  distributions.  Therefore,  com¬ 
parison  of  these  plots  to  those  of  Stringer  yield  only  marginal  evidence  of 
instabilities.  However,  of  note  are  the  depressed  temperature  ratios  between 
component  ions  and  electrons,  roughly  coinciding  to  the  region  following 
the  double  layer.  These  plots  are  also  useful  for  observing  the  velocities  of 
individual  beam  components.  At  left,  in  Figures  5.4,  we  have  the  presence 
of  a  sheath.  In  Figure  5.5,  the  top  left  panel  indicates  the  effect  of  the  po¬ 
tential  in  slowing  the  ions.  Irregularity  of  negative  velocities  in  the  middle 
left  plot  admits  to  instabilities. 

5.3.5  Electric  Field  Potential 

The  relative  charge  densities,  electric  field,  and  electric  potential  are  plotted 
in  Figure  5.6  around  upet  =  2750.  While  the  electric  field  acceleration 
remains  small  during  any  time  step,  the  cell  to  cell  variation  in  density  can 
be  quite  large.  A  smoothing  algorithm,  used  to  “best  fit”  the  data  within 
bands  of  fixed  “standard  deviation”,  resulted  in  the  continuous  curves  in 
the  density  and  electric  field  plots.  A  double  layer  is  evident  at  2  =  800 A# , 
where  there  is  a  large  potential  drop,  a  spike  in  the  electric  field,  and  the 
characteristic  back  to  back  positive  and  negative  layers  of  charge.  A  smaller 
double  layer  is  apparent  at  z  «  875Ap.  The  magnitude  of  the  maximum 
potential  drop  is  V  ~  6.5  —  7.5,  while  the  total  potential  difference  between 
particle  injection  at  z=L  and  loss  at  z=0  is  3.25.  Also  apparent  is  a  sheath 
at  2  =  0,  indicated  by  the  presence  of  the  negative  electric  field  gradient 
and  the  jump  in  potential. 

5.3.6  Wave  Spectra 

At  right  of  the  electric  field  and  potential  in  Figure  5.6  are  the  energy 
densities  in  k-space.  The  mode  structure  gives  wavelengths  most  strongly 
peaked  near  k  =  0  and  around  kX^/n  ~  .1.  These  correspond  to  longer 
scale  lengths  visible  in  the  electric  field,  l  ~  40Ao,  and  shorter  ones,  ^  ~ 
20 Ad,  apparent  in  the  potential  and  charge  density. 

Plotted  in  Figure  5.7  are  the  frequency  spectra  for  different  positions. 
These  spectra  are  measured  in  the  rest  frame  of  the  simulation  but  are 
Doppler  shifted  (uj  ±  kV)  by  the  constituents’  velocities.  The  units  of 
the  abscissa  are  tenths  of  nuipe.  The  spectra  are  most  strongly  peaked 
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Figure  5.6.  Plotted  on  this  panel  at  left  are  the  normalized  charge  density 
(top)  and  electric  field  (middle)  averaged  from  upet  =  2750  to  3000.  On  the 
bottom  is  the  electric  potential  normalized  to  the  perpendicular  injection 
temperature.  At  right  are  the  Fourier  energy  densities  El /An  (top)  and 
Pk<t>k  (bottom)  normalized  to  particle  thermal  energy  density.  The  abscissa 
for  these  plots  is  the  wavenumber  n  =  The  potential  and  spectra  are 

averaged  over  100  upeAt  centered  at  upet  =  2750. 


near  0.3  for  large  z  and  around  .06  and  its  harmonics  for  decreasing  z. 
These  correspond  to  oscillations  at  the  electron  and  ion  plasma  frequencies 
respectively. 

5.3.7  Energies 

In  Figure  5.8  we  plot  energies  normalized  to  the  injected  energies  of  the 
particles  measured  during  the  simulation.  By  u jpet  «  1250  the  plasma  has 
crossed  the  simulation  region  and  energy  is  lost  at  s=0.  Loss  rates  of  ion  and 
electron  energies  indicate  comparable  exit  currents  as  required  to  maintain 
overall  charge  neutrality  in  the  simulation.  The  collective  potential  energy 
is  a  small  fraction  of  the  total.  Midway  through  the  simulation  the  waves 
have  saturated  and  the  ion  and  electron  energies  approach  steady  state 
with  the  ion  energy  about  twice  that  of  the  electrons.  As  evident  in  the 
first  panel,  total  energy  is  conserved  to  within  much  less  than  1  percent. 


5.4  Discussion 


5.4.1  Plasma  Instabilities 


Although  a  spread  in  parallel  velocities  is  a  natural  consequence  for  mir¬ 
roring  particles  with  thermal  perpendicular  velocities,  the  observed  shape 
of  the  distributions  in  these  simulations  is  accounted  for  only  by  additional 
processes.  In  absence  of  “collisions”  the  parallel  distribution  function  for 
the  two  species  is,  using  the  notation  of  Chapter  4  (Eqn.  4.5), 
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and  the  peaks  of  the  two  species  separate  as  they  encounter  the  increasing 
quasineutral  potential.  Owing  to  the  Heaviside  function,  the  shape  of  these 
distributions  presents  a  cold  ion  beam  but  thermal  electron  distribution  to 
waves  traveling  at  speeds  greater  them  the  ion  drift  velocity  but  less  than 
that  of  the  electrons  and  ions,  oscillating  at  their  plasma  frequency,  may 
be  excited  by  electrons  drifting  through  them.  This  instability  accounts  for 
the  observed  drifting  electron  distribution  with  a  peak  coincident  to  the 
ions’,  but  with  a  “temperature”  on  the  order  of  the  electric  field  potential, 
eA< j>  ~  /iAB,  necessary  to  maintain  quasineutrality.  This  temperature 
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Figure  5.8:  Plotted  in  these  panels  is  the  energy  budget  vs.  time  ( upet ) 
for  the  simulation.  At  upper  left  is  total  particle  energy  normalized  to 
the  total  injected  energy.  This  account  is  divided  between  energies  of  ions 
and  electrons  as  displayed  in  the  two  succeeding  panels.  At  bottom  left  is 
the  normalized  potential  energy  (£,  q<f>i/Einj.)  To  its  right  is  the  energy 
contained  in  particles  lost  from  the  system. 
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spread  is  evident  both  in  the  plots  of  the  particles’  phase-space  (Fig.  5.2) 
and  their  velocity  distributions  (Fig.  5.3). 

As  the  beams  are  reflected  by  the  magnetic  field  and  electrons  are  ac¬ 
celerated  in  the  double  layer  potential,  the  constituents  become  susceptible 
to  other  two-stream  instabilities.  Electrons  in  one  beam,  excited  by  those 
in  the  opposing  beam  both  with  mean  velocities  v  >  1.3vthe,  take  part  in 
oscillations  with  u  «  uipe  and  A  ~  \d  and  acquire  a  spread  in  velocities  on 
the  order  of  the  ions’  drift  [47]. 

With  T;  <C  Tt  the  plasma  is  vulnerable  to  electron/ion  instabilities  at 
the  ion  plasma  frequency.  Figures  5.4  and  5.5  show  temperature  ratios, 
particularly  for  down-streaming  ions  on  the  high  potential  side  of  the  dou¬ 
ble  layer,  potentially  unstable  to  electron  beams  drifting  through  them  with 
velocities  greater  than  Vc  =  . 924(1. 2)v^[27],  either  combined  or  singly,  with 
two  attendant  possible  modes  of  oscillation.  Consistent  with  the  observed 
spectra,  these  oscillations  have  u>  ~  u>pi-  and  wavelengths  on  the  order  of 
Ajr>,-,  the  local  for  Te  ~  \Mv7d.  Both  ions  and  electrons  participate  in 
these  oscillations.  However,  most  of  the  excitation  energy  is  converted  to 
ion  thermal  energy  when  the  wave,  having  trapped  the  ion  peak,  breaks 
and  heats  ions  at  expense  of  the  bulk  of  electrons,  accelerated  by  the  elec¬ 
tric  potential.  Paranthetically,  the  ion-ion  wave  is  not  permitted  in  our 
simulation,  since  we  average  charge  over  a  debye  length. 


5.4.2  A  Kinetic  Approach  to  Double  Layers 


The  Langmuir  condition  requires  that  the  difference  in  electric  field  pressure 
across  the  double  layer  be  zero.  Any  difference  in  pressure  will  result  in 
movement  of  the  double  layer  until  equilibrium  is  reached.  This  equilibrium 
need  not  be  stationary. 

In  the  rest  frame  the  double  layer  momentum  is 

rx+(x,t)  , 

P  =  2^  /  Ada;  /  dvmafa(x,v,t)v 

a  •'*— (*>0  •' 


where,  as  depicted  in  Figure  5.9,  x+  and  x_  axe  the  time-dependent  spatial 
limits  of  the  double  layer  defined  by  n,  =  ne  and  A  is  its  constant  cross 
section.  The  pressure  on  its  surface  must  be 
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where  the  second  term  in  parentheses  is  the  change  in  momentum  due  to 
movement  and  changing  size  of  the  DL  while  the  first  term  is  the  pressure 
due  to  the  time- dependent  behaviour  of  the  particle  distributions. 

Ignoring  collisions  and  in  the  absence  of  instabilities,  the  Vlasov  equa¬ 
tion  3.2(a)  allows  us  to  substitute  for  the  first  term 

(^+v-Vf)  =  -a-VJ 

The  second  term  may  be  broken  into  motion  of  the  DL  with  velocity 
defined  at  its  center  where  n<  =  ne,  and  the  rate  of  change  in 
its  volume  due  to  movement  of  its  endpoints 

H  =  It(L  +  L  VDldt)  =  TtJrVDL 

The  equation  for  the  pressure  becomes 

P  =  dvmav[-  J  dxa  •  V„/  +  ( vDL  +  ^)/£t] o 
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Only  acceleration  and  reflection  of  ions  and  electrons  in  the  electric  field  of 
the  DL  contribute  to  its  recoil.  This  component  is  present  in  the  first  term 
while  the  second  term  admits  the  resulting  movement  of  the  DL. 

When  the  first  term  is  in  balance,  an  equilibrium  may  exist  with  vql  + 
^  =  0.  For  stationary  double  layers  and  symmetric  velocity  distributions 
the  second  term  vanishes  and,  in  absence  of  external  forces,  an  equilibrium 
can  only  be  established  with 


d 

53  J  dxaE  J  dvv—fa  =  qa  J  dxE(vf\t  -  J  dvf)a 


_  r  r  dE 

=  2^qa  J  dxEna  =  J  dxE—  —  0, 


where  we  have  integrated  by  parts  and  used  Poisson’s  equation  3.1(a).  This 
is  Langmuir’s  jump  condition.  In  steady  state,  for  stationary  double  layers 
and  non-symmetric  particle  distributions,  current  must  flow  to  maintain  a 
momentum  balance  between  accelerated  ions  and  electrons  [6]. 

When  one  includes  wave  particle  interactions,  an  additional  term,  §£| coii, 
contributes  to  evolution  of  particle  distributions.  The  distributions  for  an 
injected/reflected  pair  of  beams  will  no  longer  be  symmetric  since  reflected 
beams  “carry  the  memory”  of  changes  in  their  distributions  along  the  entire 
path  they  traverse.  That  is  f(x,v—)  f(x,v+).  In  this  instance  a  non¬ 

stationary  equilibrium 


53  J  dvmav(vDL  +  —  )fa\t 
or  P(vDL  +  ~)J1 
with  vDl 


=  -53/  dvmav  J  dVi^Un) 
=  RHS 

RHS  PdL/dti+ 

AP  A  P  '~ 


may  be  found.  The  physical  interpretation  is  that  the  double  layer  will 
move  or  change  in  volume  in  response  to  changes  in  particle  momentum 
[26,  p.  5,eq.  27  a,b]  until  an  equilibrium  with  vdl  +  ^  -  0  is  reached. 
Formally,  we  may  approximate  the  collision  term  as 

=  —(do  •  V„/i  +  Si  •  V„/0) 
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Figure  5.10:  Plotted  in  this  figure  is  maximum  potential  energy  vs.  injected 
ion  energy  in  terms  of  the  particle  thermal  energy.  The  analytic  expression 
for  the  expected  points  plotted  in  this  graph  is  given  by  Eqn.  4.6. 

If  the  scale  length  of  the  DL  is  far  removed  from  the  wavelength  of  the 
instability  (L  A),  we  may  ignore  the  contribution  of  the  DL  potential  to 
velocity  space  diffusion  from  instabilities,  that  is 

d/o,  _  d  q  2  E\  df0 
dt  coU  dv  m  lu  —  kv  dv  ’ 

In  our  simulation  with  ff|co«  ~  ~vf,  v  <  £ind  dx  &  A d(^)  a  rough 
estimate  of  the  equilibrium  drift  velocity  is  vdl  =  v\ d(^)  <  c„  [18],  on 
order  of  that  observed. 

5.4.3  Comparison  with  Plasma  Jet  Analysis 

In  Chapter  4  we  calculated  the  potential  due  to  a  neutral  jet  of  plasma 
with  constant  parallel  velocity  but  with  perpendicular  Maxwellian  velocity 
distributions.  Plotted  at  Figure  5.10  are  the  kinetic  energies  of  injected 
ions  vs  maximum  and  total  potentials  for  simulations  with  different  mass 
ratios.  Simulations  with  cold  plasma  injected  at  the  “ionosphere”  exhibited 
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little  difference  between  the  two  potentials  and  therefore  appear  as  single 
points. 

The  two  sets  of  lines  visible  in  the  data  axe  estimates  of  the  maximum 
potential  (solid  lines)  and  total  potential  difference  from  injection  to  exit 
(dotted  lines).  As  discussed  in  Chapter  4,  the  maximum  maintains  local 
quasineutrality  while  the  injection  voltage  assures  global  charge  conserva¬ 
tion.  Our  analysis  predicted  a  roughly  linear  dependence  between  both 
potentials  and  the  injected  ion  kinetic  energy.  Additionally,  the  total  po¬ 
tential  is  dependent  on  both  the  temperature  ratio  and  the  mass  ratios  of 
the  two  species  and  the  mirror  ratio  between  injection  and  exit.  In  the  limit 
of  large  mass  and  mirror  ratios  the  expected  slope  for  the  total  potential 
was  |  of  the  ion  kinetic  energy  for  T,  =  Te.  These  predictions  are  the  two 
lines  passing  through  the  origin. 

The  second  set  of  lines  with  non-zero  y-intercepts  are  the  least-squares 
best  fits  to  the  data.  Use  of  cold  ions  and  electrons,  effective  Te/T{  1,  in¬ 
stabilities  which  modified  particle  distributions,  and  caused  non-monotonic 
potentials  contributed  to  differences  between  the  two  sets  of  lines. 

We  compare  the  potential  expected  from  two  drifting  Maxwellians  to 
that  from  a  beam  of  ions  and  Maxwellian  electrons  by  dividing  the  potential 
into  four  regions.  In  the  first,  the  ion  density  is  approximately  that  for  a 
cold  beam  where  n,  oc  F  ~  For  quasineutrality 


nei  =  n,0  = 


e<f>i  —  fiAB 
T 


Hi  Vi  l 

"o-p- /(l  -  -jr)2 

&Q  i\0 

i  Ei _  Ii  n  _  Yl\ 

**  p  9  **  ( 1  TT  ) 

(ln|i  +  y;)/(l--^) 

■Do  0 


and,  unlike  the  case  for  drifting  Maxwellians,  the  potential  in  region  I  is 
only  weakly  dependent  on  the  ion  kinetic  energy. 

The  boundaries  of  the  second  region  are  defined  by  points  just  before 
and  after  the  double  layer.  Noting  that  the  ion  density  may  increase  from 
its  injected  value  before  an  inflection  point  at  F(1.5)  =  .43  to  a  maximum 
defined  by  F(. 92)  =  .54  [1,  p.298],  a  rough  estimate  of  the  magnitude  of 
the  DL  is  obtained  by  setting  the  argument  of  F  to  1. 


Ki  -  VIt 
V/ 


1 
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Vi  i  «  Ki-Yi 


and  the  maximum  potential  is  Vm  =  V\  +  Vu  w  Kj  +  In  The  value  of 
the  potential  drop  is  directly  dependent  on  the  ion  kinetic  energy  at  the 
DL’s  base.  Indeed,  since  the  DL  must  begin  before  the  inflection  point, 
<  2.25  or  Ki  >  Vu  >  Kt  -  2.25 Yi  and  F/  <  Ki /2.25.  The  fact 
that  the  potential  is  less  than  the  injected  ion  kinetic  energy  relates  to  the 
investment  of  ion  energy  into  thermal  energy  of  the  electrons. 

One  key  to  double  layer  existence  is  the  slope  of  the  electron  distribution 
function  vs  that  of  the  ions.  Again,  using  n,  oc  F, 

K0-V2  j  <  2{M)X, 

riu  Yu  Yi 

Similarly,  quasineutrality  demands 


rieii/nei  <  2(.54)A7 

exp1'"  <  1.08X; 

Vu  <\nXj 

<  1 1  Kl 

-  2  n  T7 

This  condition  is  equivalent  to  a  Bohm  criteria  where  the  electron  slope  at 
the  double  layer’s  start  is  required  to  be  greater  than  the  ions’  (^  >  ) 

but  less  than  that  at  the  inflection  point  Satisfaction 

of  a  similar  relationship  at  the  high  potential  end  of  the  double  layer  is 
guaranteed  by  the  shape  of  the  ion  distribution  function  [35]. 

Because  electrons  gain  and  ions  lose  a  substantial  fraction  of  the  ion 
kinetic  energy  at  the  double  layer,  in  the  third  region  a  decreasing  poten¬ 
tial  may  be  required  to  maintain  quasineutrality.  Generally,  however,  the 
potential  is  relatively  flat  (V//j  =  0)  and  we  approximate  both  ions  and 
electrons  as  beams. 

In  the  fourth  region,  at  the  boundary  of  the  simulation,  a  sheath  must 
exist  to  balance  the  electron  and  ion  current  losses.  For  cold  beams  a  rough 
estimate  of  the  needed  potential  is 


Viv  = 


Ku  -  kn 


Vi  -  (Vii  +  h) 

2 

Vi  -  (K,  -  Yt)  -  kj  ,J0  +  */ 
- n - =  -( — o - Yn 


so  the  total  potential  for  current  equality  is 


Vi  +  Vu  +  Vm  +  VJV  =  (\n-j-  +  Yi)  +  (KI-YI)  +  0-(KIl2-YI  +  kI/2) 

«  In  +  Yi  +  Ki/2  -  hi  12 

-DO 

While  the  intercepts  of  the  two  potentials  are  sensitive  to  double  layer 
location  and  the  actual  distributions  of  ions  and  electrons,  the  relative 
slopes  between  the  two  potentials  observed  in  our  simulations  agree  well 
with  this  heuristic  approach. 
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Chapter  6 


2D  Simulation 


Results  for  two-dimensional  runs  are  presented  for  comparison  to  the  one¬ 
dimensional  results  of  the  previous  chapter.  Although  the  physical  model 
is  the  same  as  for  ID,  the  two  dimensional  simulation  requires  extending 
the  Poisson  solver  beyond  simple  integration  used  in  the  previous  chapter. 
We  discuss  these  aspects  more  fully  in  Appendices  B  and  C.  Similarly,  as 
discussed  in  Appendix  A,  refinement  of  the  particle  pusher  is  required  to 
permit  following  greater  numbers  of  particles  contained  in  our  two  dimen¬ 
sional  “volume”  for  a  reasonable  period  of  time. 

6.1  Computational  Model 

6.1.1  Electric  Field 

Electric  Field  Grid 

The  model  and  boundary  conditions,  pictured  in  the  accompanying  sketch, 
are  similar  to  those  in  z  for  ID,  but  now  the  electric  potential  also  varies 
radially  in  r.  Now,  to  collect  charge  and  compute  the  potential,  a  two- 
dimensional  mesh  in  r  and  2  overlays  the  simulation  region.  The  cell  size 
for  both  directions  is  \q,  as  in  the  ID  case. 

Charge  Sharing/Force  Interpolation 

The  weighting  scheme  chosen  in  ID  may  be  interpreted  as  a  dipole  expan¬ 
sion  of  the  charge  density  about  the  midpoint  between  cells  [29].  In  two 
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Figure  6.1:  The  two  dimensional  model.  The  coordinates  at  top  left  are 
z,r=0.  At  bottom  right  they  Eire  z=L,  r=R.  Boundary  conditions  in  z  are 
those  for  ID.  That  at  r=Q  are  appropriate  for  continuous  charge  distribu¬ 
tions.  Although  we  show  <f>(R)  =  0,  the  boundary  condition  at  r=R  is 
'lexible. 


dimensions  a  similar  expansion  including  a  quadrupole  term  yields  a  scheme 
known  as  area  weighting  [5,  p.  244] .  This  approach  is  roughly  equivalent 
to  the  scheme  used  here. 

In  ID  charge  sharing/force  interpolation  derived  from  treating  particles 
as  plates,  one-cell  thick  of  uniform  charge  density.  In  two-dimensions  (r-z) 
we  may  use  the  same  approach  in  the  z  direction.  In  r,  however,  we  treat 
the  particle  as  a  one  cell  thick  ring  which  varies  radially  to  maintain  total 
charge,  q  —  ir p(r2,t  —  r2  JAz.  Therefore, 

l+2  1  2 


P  = 


2irr,  AzAr 
,—.,2 _ 


■+i 


<1 
^  2 


where  r,-+i  and  r,_i  are  the  outer  and  inner  bounds  of  the  particle. 
Using  Gauss’  Law,  the  electric  fields  at  a  radial  location,  r,  are 


rET 


r+[6Er(r+,z+)  +  (1  -  6)Er(r+,z.)}(^ - 

r+  —  rt 

2  2 

+r_[«£r(r.,2+)  +  (l-«E.(r.,z.)](^^T) 

r+  —  r  i 


Ez  =  ^[^(r+,z+)( 


r2  —  r2 


r2  _  r2 


)  +  Ez(r-.,z+)(: 


0] 


+  (I  -6){E2(r+,z_)( 


r2  —  r2 
r2  —  r2_ 


ri  —  r2 


)  +  ^(r_,z_)(-f~r)], 
r+  —  ri 


where  the  pluses(-f)  and  minuses(-)  indicate  the  upper  and  lower  cells  into 
which  the  particle  extends. 

The  radial  and  axial  electric  forces  are  then  calculated  as 
Fr  =  2  n  rdrpET  =  q+ET+  +  q~Er_, 

—  =  EZ  =  6Ez++(\-6)Ez_ 

q 

One  must  use  care  in  calculating  the  force  due  to  the  magnetic  field. 
This  force  should  act  at  the  charge  center  of  the  particle  [8]  which  is  located 
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at 


rc 


(ri  +  £ )/\/2  r,  <  \ 


Poisson  Solver 

Computation  of  the  electric  field  is  obtained  by  implementing  Buneman’s 
algorithm  for  mixed  boundary  conditions  (Appendix  C),  solving  for  the 
potential  and  ultimately  the  electric  fields. 

Boundary  Conditions 

The  boundary  conditions  at  z  =  0  and  L  are  those  of  the  preceding  chap¬ 
ter.  The  boundary  condition  at  r=0  attests  to  the  continuous  particle 
distributions-no  line  charge.  For  high  order  radial  modes  the  boundary 
condition  at  r=R  is  flexible.  In  two  dimensions  we  may  express  the  solu¬ 
tions  to  Poisson’s  equation  3.1(a)  in  terms  of  Bessel  functions  (Appendix 
B )-(f>(R)  =  0  or  Er(R )  =  0  just  shifts  the  solution  by  a  node.  Since  we 
model  a  wide  enough  radial  region  to  expect  the  presence  of  several  modes, 
this  choice  is  not  critical. 

At  r  =  R  and  z  =  0  we  allow  the  charge  in  particles  “hanging  over  the 
edge”  to  be  lost.  At  z  =  L  we  take  advantage  of  inversion  symmetry  and 
mirror  particles  as  for  ID. 

6.1.2  Particles 

Particle  Pusher 

For  the  ions  we  retain  the  full  Boris  pusher  of  the  previous  simulations. 
For  the  electrons,  however,  we  adopt  a  guiding  center  pusher,  as  described 
below.  In  particle  pushing  schemes  including  gyromotion,  length  of  the  time 
step  is  constrained  by  the  gyroperiod  (Appendix  A.)  In  a  simulation  with  a 
magnetic  field  gradient,  such  as  a  dipole’s,  this  limits  the  magnitude  of  the 
magnetic  field  at  exit,  and,  when  the  field  at  injection  assumes  particular 
values,  ultimately  the  mirror  ratio.  We  desire  to  adopt  a  scheme  which 
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ameliorates  this  condition.  Averaging  particle  motion  over  Larmor  orbits 
results  in  the  set  [32] 

mdv  ||  pdB  m_  (dt\  dtx  . 

7  a  =  "  7  aT  +  ' ( aT  + +  U£ ' Ve,) 

i  ei  ,  ->  uc- ♦  me.  de'i  duE„ 

R1  =  -x{-cE  +  -VB  +  — +  -^-1} 

The  desirability  of  time-centered  schemes  is  discussed  in  Appendix  A. 
In  general,  this  set  of  equations  is  not  amenable  to  time- centering.  However, 
in  cylindrical  coordinates,  because  both  E  and  stationary  B  lie  in  the  /f,  z 
plane, 

dv H  eE\\  ii  dB  .  . 

-Jl  =  —II  _  •  uE  ■  V  6 

dr  m  m  os 


v±  =  r<j><f>  =  Vd  + 


q  <  E  >  xb 


where  ^  is  entirely  in  the  ^ 

direction  for  low  oscillation  frequencies  (w  <C  ft)  and  small  radial  scale 
lengths  (p  <C  f2e).  With  axial  symmetry  this  perpendicular  component  can 
be  ignored  in  moving  particles. 

Subtracting  v±  from 

dx  *  a 

— -  =  rr  +  r<t><t>  +  if 


leaves 

{T||  =  rr  -f  ii 

and  differentiating  the  parallel  velocity  leads  to  a  time-centered  expression 
in  the  r-z  plane, 

*4 -rf  +  zz  -  ^h  +  v2(^r  +  ?^z) 
dt  ~rr  +  zz  ~  dtb  +  V"{dsr+  dsZ) 

— r  +  ^  — — 2  i  =  -(“  3uj|(l  +  cos2 $)be(b  x  ^)/r/(l  +  3cos2  0)]o 

Strictly,  the  guiding  center  approximation  is  not  applicable  within  the 
double  layer  itself.  However,  as  demonstrated  in  one  dimension,  for  suffi¬ 
ciently  narrow  double  layers,  the  pusher  well  models  the  bulk  of  the  simu¬ 
lation  volume. 
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Particle  Injection 

The  particles  are  distributed  in  r  and  injected  uniformly  at  the  injection 
boundary.  As  for  our  ID  simulation  the  full  dipole  magnetic  field  is  super¬ 
imposed  on  our  mesh.  However,  for  the  two  dimensional  model  the  dipole 
is  mirrored  at  the  injection  boundary  ( z  =  L),  and  the  symmetric  field 
calculated  accordingly.  This  technique  allows  injection  entirely  in  the  z  di¬ 
rection  while  distributing  the  particles  in  r.  The  radial  injection  is  designed 
to  extend  many  ion  Larmor  radii  while  keeping  plasma  from  being  lost  out 
the  sides. 


6.2  Two  Dimensional  Results 

6.2.1  Simulation  Parameters 

The  attached  plots  are  for  a  mass  ratio  of  25,  a  mirror  ratio  7  ~  11, 
simulation  length  of  L  =  128Ao  and  width  R  =  64A d-  Temperatures  are 
Ti  =  2Te  =  lOOOfceF.  Both  ions  and  electrons  are  injected  uniformly  with 
vdrift  ~  3\Z2/bvthe  at  the  boundary  z  =  L  for  r  <  32 A^.  These  are  the 
same  values  as  for  the  previous  simulations  and  therefore  have  the  same 
relative  scale  lengths  and  times. 

6.2.2  Chronology 

As  pictured  in  Figure  6.2,  the  electrons  are  reflected  rapidly  by  the  mirror 
field.  Initially  the  electrostatic  potential  is  ~  fcTe,  as  the  ions  attempt  to 
neutralize  the  front-going  electrons.  As  for  the  one- dimensional  case,  they 
rapidly  attain  an  effective  parallel  temperature  ~  kTe.  Particles  with  large 
drift  velocity  may  result  from  E  x  B  drifts. 

6.2.3  Velocity  Distributions 

Underlying  velocity  space  plots  for  both  ions  and  electrons  are  oscillations 
of  amplitude  <  kTe.  As  in  the  one  dimensional  case,  ions  fill  velocity  space. 
The  profile  of  the  velocity  distributions,  Figure  6.3,  moderately  into  the 
injection  region  shows  ion  velocities  embedded  in  a  two-  peaked  electron 
velocity  distribution.  The  relative  drifts  between  the  electrons  not  only 
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Figure  6.2:  Plotted  in  this  chapter  are  run  results  for  ui^t  —  200  (top)  and 
wpet  =  250  (bottom).  In  this  figure  the  left  plots  are  normalized  ion  (top) 
and  electron  (bottom)  velocities,  the  middle  set  is  average  velocities,  and 
at  right  are  fractions  of  trapped  (reflected)  particles  versus  axial  position. 
Ion  velocities  in  the  left  panels  are  normalized  to  half  the  electron  thermal 
velocities  (upeAt  =  .5)  to  which  the  remainder  are  normalized. 
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allow  the  ions  to  interact  with  both  distributions  but  may  be  sufficient  for 
an  electron-electron  two-stream  instability. 

While  use  of  the  guiding  center  pusher  eliminates  opportunity  to  ob¬ 
serve  electron  cyclotron  waves  (partially  responsible  for  diffuse  aurorae), 
one  would  expect  electrostatic  ion  cyclotron  modes  (fcx  ^  &||)  to  be  present. 
Even  though  the  simulation  proceeds  for  many  DA t,  the  influence  of  these 
waves  cannot  be  distinguished  from  these  plots.  Overall,  the  gross  results 
compare  well  with  the  ion  acoustic  scenario  presented  in  the  previous  sec¬ 
tion. 

6.2.4  Electric  Field  Potential 

The  large  scale  potential  is  of  the  order  to  maintain  charge  neutrality.  At 
later  times,  uptt  ~  200,  a  sharp  potential  drop  is  readily  apparent  from 
s  ~  75.  The  magnitude  is  ~  7 kTe.  However,  its  movement  is  not  apparent, 
or  it  actually  may  be  receding.  Reminiscent  of  the  V-shaped  potentials 
encountered  in  nature  [51],  the  potential  contour  plots  clearly  evidence  a 
gross  radial  scale  size  of  the  jet  ~  32 Ad  (as  broadened  by  the  ion  gyroradius) 
necking  down  to  about  1/3  that  value  at  the  left  boundary,  corresponding 
with  a  similar  increase  in  B.  The  local  scale  size  for  the  potential  is  far 
more  granular. 

Also  plotted  here,  the  densities  show  a  marked  increase  as  the  particles 
slow  down  in  the  effective  potential  and  then  slowly  decrease  to  the  left 
boundary  as  particles  are  reflected. 


6.3  Discussion 

Quantitatively  and  qualitatively  the  two-dimensional  simulation  results  agree 
well  with  the  one- dimensioned  results.  This  is  evidenced  particularly  in  the 
velocity  space  diagnostics.  Although  we  may  expect  break-down  of  the 
guiding  center  approximation  for  /  Et  dt  too  great  and  from  ignoring  com¬ 
ponents  of  particle  movement,  gross  simulation  energy  is  well  conserved,  as 
shown  in  Figure  6.5.  The  question  of  stationarity  in  2D  is  open.  While, 
the  solutions  may  be  altered  by  the  additional  source  term  (f;f)  in  the 
Poissons  equation  from  its  ID  form  [9],  the  additional  degree  of  freedom 
for  particles  and  waves  need  further  exploration. 
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Figure  6.3:  Plotted  are  the  velocity  distributions  for  the  electrons  and  ions 
at  eight  points  from  2=0  (top  left to  z  =  L  (bottom  center). 


Figure  6.4:  Plotted  in  this  figure  are  the  particle  densities  (at  left),  relative 
charge  number  (top  center),  and  scale  length  (bottom  center).  The  right  set 
of  plots  portray  the  potential  along  the  axis  and  radially  by  axial  position. 
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Figure  6.5:  These  plots  show  the  same  energies  for  the  2D  case  as  for  the 
similar  plots  in  ID  (Figure  5.8). 
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Chapter  7 
Conclusion 


Our  dissertation  goal  was  to  investigate  the  detailed  spatial  and  temporal 
behaviour  of  potentials  in  the  scenario  of  Serizawa  and  Sato  [42].  We  ob¬ 
tained  analytic  confirmation  of  their  semi-empirical  result  that  global  large 
scale  potentials  may  exist  for  plasmas  injected  into  a  dipole  magnetic  field. 
The  potential  is  proportional  to  the  kinetic  energy  of  the  injected  ions,  as 
given  by  equation  (4.6).  However,  our  analyses  in  Section  4.3  further  sug¬ 
gest  that  local  double  layer  solutions  are  to  be  expected  in  the  parameter 
regimes  considered.  These  expectations  were  verified  by  particle  simulation 
in  one  and  two  dimensions  and  the  results  are  described  in  Chapters  5  and 
6. 

We  were  led  to  our  analytic  results  in  search  of  a  simple  computational 
model.  We  therefore  investigated  analytic  properties  of  increasingly  more 
realistic  particle  distributions.  In  a  cold  drifting  plasma  the  current  and 
density  are  locked  together  by  equation  (4.4)  so  that  the  global  and  local 
solutions  are  inextricably  linked.  Thus,  in  Section  4.1  we  demonstrated 
that,  although  a  large  scale  potential,  equation  (4.2),  exists,  it  neither  per¬ 
mits  local  double  layer  solutions  nor  guarantees  acceleration  of  incident 
particles. 

With  perpendicular  and  parallel  temperatures,  global  potentials,  Eqn. 
(4.6),  retain  the  same  functional  dependence  on  ion  kinetic  energy  (~  if,) 
as  for  the  cold  case,  and  maps  of  their  relative  charge  densities  exhibit  re¬ 
markably  constant  properties.  However,  some  distributions  lead  to  local 
double  layer  solutions,  others  do  not.  In  sections  4.2  and  4.3  we  demon¬ 
strated  that  the  local  double  layer  solutions  are  provided  opportunity  by 
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the  Maxwellian  plasmas  and  cold  (parallel)  ions.  Local  analyses  indicated 
this  in  section  4.3.2  and  Figure  4.6. 

The  analytic  results  can  be  couched  in  terms  of  global  and  local  require¬ 
ments  for  charge  neutrality.  Global  neutrality  is  guaranteed  by  equal  exit 
currents  into  the  loss  cone  while  demands  of  local  neutrality  may  lead  to 
double  layers.  The  conditions  for  a  double  layer  are  indicated  by  multi¬ 
valued  solutions  to  the  quasi-neutral  equation  (3.6).  These  solutions  are 
present  only  when  the  slopes  of  underlying  velocity  distributions  permit 
(as  allowed  by  Bohm’s  criteria  discussed  in  Section  3.4.1)  and  the  mo¬ 
mentum  balance  of  particles  rebounding  or  accelerating  in  the  potential 
is  maintained  (Langmuir’s  condition.)  Equivalently,  Langmuir’s  condition, 
equation  (3.9),  may  be  interpreted  graphically  as  areas  of  equal  charge 
densities  integrated  over  the  double  layer  potential  (Section  5.4.2.) 

Indeed  Bohm  criteria  and  Langmuir’s  condition  are  just  as  operative  as 
ever.  All  the  theory  of  Chapter  3  is  consistent.  In  section  3.4  we  showed 
that  the  Bohm  criteria  were  equivalent  to  conditions  on  the  charge  density. 
In  section  3.6  we  related  these  conditions  to  treatment  of  currentless  dou¬ 
ble  layers.  However,  the  Bohm  criteria  are  necessary  but  not  sufficient  to 
guarantee  a  solution.  While  the  frame  of  the  double  layer  is  important,  the 
real  requirement  is  on  particle  distributions,  not  average  quantities  such  as 
injection  velocity. 

To  explore  the  temporal  behaviour  of  our  analytic  results,  in  Chapters 
5  and  6  we  simulated  a  drifting  plasma  incident  upon  a  dipole  magnetic 
field  in  one  and  two  dimensions.  While  previous  double  layer  simulations 
modeled  constant  magnetic  fields  and  fixed  potentials  [43],  to  our  knowledge 
this  is  the  first  simulation  of  DL  formation  combining  particle  injection, 
floating  potentials,  and  realistic  magnetic  fields.  In  the  one  dimensional 
simulations,  a  double  layer  of  several  kT  with  scale  length  /  ~  A©  <gC  R,  the 
magnetic  field  scale  length,  is  clearly  apparent  in  Figure  5.6. 

As  discussed  in  section  5.4.3  and  shown  in  Figure  5.10,  our  simulation 
results  compare  well  with  analytic  predictions  [42],  and  correspond  qualita¬ 
tively  with  laboratory  experiments  [44],  theory,  and  simulations  examining 
the  role  of  plasma  jets  into  the  earth’s  polar  regions  [42]  and  agree  in  mag¬ 
nitude  with  their  predictions.  In  contrast,  previous  theory  [52]  and  simu¬ 
lations  [14]  predicted  an  electric  field  E  a  jL  (/  ^  j?.)  These,  however, 
assumed  self-consistent  distributions  of  particles  in  quasineutral  electro¬ 
static  fields,  thereby  ignoring  those  portions  of  solution  space  reserved  for 
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double  layers.  Our  results  bear  out  Stern’s  more  general  approach[45]  of 
treating  the  double  layer  as  an  extension  of  quasineutrality  but  on  a  coarser 
(than  Ad)  scale. 

Our  simulations  have  proven  valuable  in  considering  access  to  double 
layer  solutions  which  match  the  required  global  potential [42].  While  the¬ 
ory  shows  that  double  layers  may  form  with  self-consistent  particle  distri¬ 
butions,  consideration  of  overall  charge  neutrality  demands  that  the  total 
potential  drop  across  the  model  boundaries  be  constrained  to  a  particular 
value.  Local  and  global  solutions  need  not  match.  Adjudication  of  any 
difference  between  these  opposing  potentials  obeys  the  spatially  depen¬ 
dent  solution  of  Poisson’s  Equation  coupled  to  the  particles’  distributions 
throughout  the  simulation  region.  While  the  global  condition  depends  on 
the  presence  of  particles  in  the  loss  cone,  the  double  layer  solution  ulti¬ 
mately  depends  on  local  conditions.  In  our  simulations  differences  between 
the  boundary  solution  required  for  global  charge  neutrality  and  the  local 
solution  for  the  potential  resulted  in  a  sheath  at  the  exit  boundary. 

Potential  and  particle  distributions  must  be  self-consistent  in  time. 
Perkins  and  Sun  [35]  analyzed  the  instability  criteria  for  distributions  such 
as  ours.  As  expected,  both  ion  and  electron  beams  are  subject  to  streaming 
instabilities.  The  instabilities,  evident  in  temperature  and  velocity  space 
plots  of  Figures  5. 1-5.5,  are  consistent  with  the  ion  acoustic  instability  (with 
both  electron  beams)  in  one  dimension.  At  first,  the  electrons  rapidly  ther- 
malize  while  the  ions  remain  cold.  The  waves  scatter  the  ions  in  velocity 
space.  Coincident  with  this  scattering,  the  double  layer  moves  to  the  low 
field  side.  The  double  layer  approaches  and  rests  near  the  boundary  (z=L). 

In  consonance  with  [35]  and  suggested  by  [6],  the  frame  of  the  double 
layer  is  important.  When  stability  is  reached  the  final  distributions  must 
be  consistent  with  the  final  potential,  and  instabilites  affect  both.  The  role 
of  instabilities  is  to  alter  particle  distributions,  and  this  requires  the  dou¬ 
ble  layer  to  move.  In  section  5.4  we  explained  the  movement  in  terms  of 
momentum  conservation,  or  Langmuir’s  condition,  versus  Block’s  descrip¬ 
tion  in  electrical  terms.  Momentum  conservation  demands  that  the  double 
layer  adjust  itself  to  where  a  local  solution  exists.  The  momentum  lost  is 
attributed  to  waves,  in  analogy  with  quasilinear  theory  [26].  Unlike  ours 
Stenzel  et  al’  double  layers  were  stationary [44].  However,  they  attributed 
this  property  to  a  species  of  trapped  electrons  characteristic  of  their  exper¬ 
iment. 
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The  theory  we  have  reviewed  is  essentially  ID  but  requires  considera¬ 
tions  of  a  second  dimension.  The  two  dimensional  code  has  been  sucess- 
fully  employed,  but  2D  results  remain  quite  preliminary.  It  shares  the  same 
global  behaviour  but  not  all  details  of  the  local  behaviour  of  the  ID  results. 
In  particular  the  2D  results  lead  to  smaller  magnitude  double  layers,  and 
their  movement  is  not  apparent.  Because  our  2D  diagnostics  were  not  as 
extensive  as  in  one  dimension,  we  are  unable  at  this  time  to  verify  the  role 
of  electrostatic  ion  cyclotron  waves.  However,  the  code  should  see  future 
utility  in  more  fully  observing  the  effect  of  the  extra  degree  of  freedom.  For 
example,  in  two  dimensions  it  may  be  useful  to  model  off-axis  behaviour 
for  our  scenario  or  to  model  other  scenarios  such  as  that  of  Perkins  and 
Sun.  We  discuss  below. 

Our  B  field  provided  an  adjustable  parameter  called  for  by  [35].  How¬ 
ever,  we  modeled  a  sudden  entrance  of  plasma  into  a  region  to  simulate 
behaviour  during  a  substorm.  In  our  scenario  we  observe  a  rapid  injection 
with  few  bounce  times.  The  distributions  are  dominated  by  injected  parti¬ 
cles.  Other  situations  may  allow  many  bounce  times,  such  as  when  particles 
are  injected  into  a  magnetic  mirror  [15].  In  this  model  both  ions  and  elec¬ 
trons  are  Maxwellian,  because  the  ions  have  time  to  explore  adjacent  cells 
where  collisions  predominate.  In  our  model  instabilities  we  observe  play  a 
similar  role.  Given  more  than  a  few  bounce  times  or  increased  collision  fre¬ 
quency  and  sufficient  confinement  time,  we  expect  to  similarly  populate  our 
model,  in  which  case  an  equilibrium  with  stationary  double  layers  should 
be  reached  [35],  [15].  Under  these  conditions  a  particle  undergoing  many 
bounces  will  lose  memory  of  its  injection  velocity.  Ultimately,  however,  the 
Langmuir’s  condition  must  be  satisfied,  requiring  the  double  layer  to  adjust 
to  some  stationary  equilibrium. 

Longer  simulation  runs  may  be  needed  to  explore  the  full  behaviour  of 
the  double  layers  such  as  these.  In  that  regard,  usefulness  of  a  guiding 
center  pusher,  developed  in  Appendix  A,  was  demonstrated.  Even  in  two 
dimensions  it  requires  no  corrector  for  stability.  However,  a  more  rigorous 
implementation  of  the  guiding  center  approximation  would  be  most  useful 
for  considering  perpendicular  phenomena  in  the  E,B  plane. 

Simulations  depending  on  Fourier  transforms  to  solve  for  the  potential 
[17]  may  be  plagued  by  errors  in  the  field  calculation.  In  our  ID  simulations 
these  errors  were  traced  to  aliasing  in  a  bounded  plasma.  In  Appendix  B 
a  correction  for  aliasing  has  been  derived  to  ameliorate  its  effects,  and 
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this  correction  was  sucessfully  employed  in  one  dimension.  Our  approach 
suggests  an  “optimal”  particle  shape  based  on  two-body  correlations  to 
account  for  collisions.  In  Appendix  C  we  developed  Buneman’s  algorithm 
for  the  mixed  boundary  conditions  of  our  bounded  2D  simulation. 

In  our  simulations  we  used  a  simple  injection  scheme  to  avoid  “energy 
instabilities.”  This  is  analgous  to  and  suffers  from  the  save  problems  as 
“quiet  starts”  [5]  to  distribute  particles  in  angles.  It  proved  adequate  for 
the  ID  results.  However,  random  injection  with  varying  velocities  (and 
densities)  of  particles  would  be  more  realistic  and  still  avoid  some  of  the 
cautions  for  steady  injection. 

Ideally,  if  charge  neutrality  is  to  be  observed  throughout  the  simulation 
region,  f  ds( =  — J )  =  0  at  the  boundaries.  Writing  this  law  in  terms  of 
the  total  charge  entering  a  flux  tube  originating  at  the  earth,  we  must  have 
/  J||  •  dAenj  =  —  j  J±  •  dAaidea •  The  influx  of  current  into  the  magnetic  field 
flux  tubes  must  be  balanced  by  that  across  the  magnetic  field.  If  there  is 
to  be  a  net  current  out  of  the  simulation,  then,  to  be  totally  self-consistent, 
one  needs  to  treat  it  as  part  of  a  circuit  as  did  Sato  and  Okuda  [37].  For 
discussion  see  Reference  [43]. 

Finally,  the  role  of  background  plasma  in,  for  example,  adjudicating 
the  potential  difference  between  local  and  global  conditions  has  yet  to  be 
fully  explored.  Its  role  was  discussed  to  some  extent  in  section  5.2.3.  Such 
techniques  are  worthy  of  future  consideration. 
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Appendix  A 
Particle  Pushers 


A.l  The  “Leap  Frog”  Scheme 


A  non-relativistic  charged  particle  in  external  electric  and  magnetic  fields 
obeys  the  Lorentz  force  law 

F  =  mf  =  ,(E+2^)  (3.3) 


[24,  p.572].  To  advance  particles  with  discrete  time  steps  (A <),  we  seek 
a  difference  equation  analog  to  (3.3)  together  with  the  similar  expression 
for  the  velocity,  &  =  v,  which  accurately  models  particle  trajectories  and 
is  also  stable  to  small  deviations  (due  to  round-off  error,  etc.)  from  the 
correct  solution. 

A  MacLaurin  expansion  of  (3.3)  gives  to  lowest  order  in  At 


v  i  —  v_  1 
2  2 


— IE>  +  C’+°~h)  X  ^  ]  +  «(Ai2) 


(A.l) 


This  equation  along  with  its  position  counterpart, 


£i  ~  Jo 

At 


is  known  as  the“Leap-Frog”  scheme.  Its  name  derives  from  alternate  deter¬ 
mination  of  velocities  and  functions  of  position,  as  pictured  in  Fig.  A.l.  We 
shall  note  that  this  time-cenb  red  scheme  provides  accuracy  and  favorable 
stability  properties. 
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Figure  A.l:  Pictured  here  is  the  Leapfrog  scheme  in  which  velocities  and 
functions  of  position  are  determined  at  interleaved  time  steps.  (After  [5, 
Fig.  2-4a,  p.13].) 


When  E  and  B  fields  are  constant,  we  may  obtain  the  solution  to  (3.3) 
by  substituting  v  =  u  +  we  where  we  =  c^§r-,  leaving 

du  q  r(u  XB)  ,  5l 

Tt  =  +  £"]’ 

where  E\\  =  E  ■  B/Bb.  Setting  Q  =  £ B/c  and  operating  on  both  sides  with 
xO  gives  =  —  Q2u±  where  u±  =  u  —  if||  and  uj|  =  u  •  B/Bb.  Using  the 
initial  conditions,  u2 o  =  U2(0)  =  0  and  U30  =  ^3(0)  =  u_lo?  the  continuous 
solution  for  the  perpendicular  velocity  is  u±  =  ux0[cosQt3  —  sinfl<2]. 

The  solution  to  (A.l)  is  obtained  in  a  similar  fashion.  To  lowest  order, 


At 


1 

2 


ui  4-  u_  1  _  _ 

— ~ — —  —  x  fio  +  E]]0 


Multiplying  each  side  by  At  and  adding  2 u_i  yield 


ul  +  u_i 
2  2 


Ul 

2 

where  a  = 


Ui  +  u_£ 
2 


—  [ui  +  u 
1  2 


x  VloAt  -I-  2u_i  -f*  E^oAt 
ijxQ;2  +  u_i  x  a  +  2u_i  +  E\\0At 


l  -  a2  2 a 

r— — ~u±_l  +u_ i  x  — — -a  +  «||_i  +  E\\0, 
1  +  a2  *  s  1  +  a2  11  2  11 
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3  =  1  x2 


Figure  A. 2:  This  is  the  magnetic  coordinate  system. 


Separating  this  (vector)  equation  into  its  three  scalar  components  ac¬ 
cording  to  the  coordinate  system  defined  in  Figure  A. 2, 

=  u^i+E^At 
u2i  =  /l«2— i  +/2«3-i 

u3*  =  hu3-%  ~  hu2-%i 

where  f\  =  and  /2  =  y^yj .  The  last  two  equations  may  be  most 
conveniently  expressed  as  u±i  =R  where  R  is  the  tensor 


/i  h 

fi  fi 


We  note  that  the  determinant  |  R  |  =  ff  +  f%  =  (y+p-)2  +  (y+p-)2  =  1  so 

4-*  — * 

the  operator  R  is  a  pure  rotation  in  the  plane  perpendicular  to  B.  This 
means  that  the  magnitude  of  the  velocity  component  perpendicular  to  B 
is  conserved. 

Because  our  particle  maintains  its  perpendicular  velocity  in  constant 
E  and  uniform  B  fields,  its  guiding  center  moves  at  the  correct  velocity, 


110 


vgc  =  t>||ei  +  we,  and  along  its  correct  path.  However,  corrections  must 
be  introduced  to  maintain  the  trajectory  for  the  particle.  Using  half-angle 
formulas  for  the  sine  and  cosine  [41],  after  one  time  step  the  continuous 

phase  angle  is  fiAt  =  tan~'(-  )•  However,  since 


u  21 

«31 


2a 


,  .  2U° 
1  +  a2 

1  —  a2 
7— — jtio, 
1  -(-  a2 


the  discrete  phase  angle  is  tan-1(— j^-)  =  While  equal  in  the 

limit  At  — >  0,  only  a  =  tana  maintains  the  same  phase  for  all  At. 

Similarly,  for  «x  to  represent  the  correct  particle  orbit,  we  must  equate 
the  difference  equation 

Pi  ~  Po 


to 


At 

dp 

dt 


=  ui 
2 


=  u 


where  here  p  is  the  Larmor  radius.  Using  the  known  solution  for  the  correct 
particle  position,  p  =  p(0)(cosftf2  +  sinQf3),  the  discrete  solution  is 

_Y  »  v  -Yr>\  —  a  ±  n  r  .  At  *  ft  At «. 

p(At)  —  p[  0)  =  u  i  At  =  fipo[—  sm  ——2  +  cos - 3]At 


giving  the  particle  position  after  one  time  step  as 

...  .  .fiAt.  _  . 

p2(A<)  =  po-  /30sin(— )StAt 

/  a  ,fiAf  . 

po\At)  =  po  cos(  — — — )fi  At 
£* 

with  Larmor  radius 


Pi  =  \jph  +  Pl\  =  Pofl  +  4a2(l  —  ———)]* 

a 

This  agrees  with  the  continuous  solution  only  for  sin  a  =  a. 

Recognizing  that,  in  general,  the  solution  to  the  difference  equation 
does  not  match  the  physical  solution,  we  may  introduce  coefficients,  A(At), 
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into  the  difference  equations  to  obtain  identical  continuous  and  discrete 
trajectories  for  all  A t  while  the  difference  and  continuous  equations  are  the 
same  in  the  limit  of  small  time  steps,  limA<_oA  =  1.  Defining  \p  = 
and  An  =  we  have 

A.LUii  + 

P2  -*-2  C2  1)2 

-  [A,  +  A„o(B„x  +  <  **— xB°)]  ( A.2) 

m  "  c 

where  <  v0  >=  ~z.  A  particularly  efficient  implementation  of  this 

scheme  breaks  the  equation  for  the  velocity  into  an  acceleration  by  the 
electric  field  and  a  pure  rotation  by  the  magnetic  field,  and  is  known  as  the 
Boris  algorithm  [8]. 

Accuracy 

— ♦  — * 

This  scheme  was  derived  for  constant  E  and  B  and  is  the  correct  solution 
for  the  first  RHS  term  in  (A.l).  Setting  x  =  x0  -f  Aar  +  e,  we  obtain  the 
error  after  one  time  step. 

X.l  -  Xq  — 

2 

Ax  = 

t  — 


L 


±iA  t 


dtv 


o 

VqA  t 


J 


vdt  —  Ax 


& 

24 


V1 


X 1  Xq 

At 
—  v 


At 


)  = 


Stability 

_>  —* 

For  simplicity  consider  (A.2)  only  for  constant  B  and  one- dimensional  E  = 

J5J||(s)6.  Setting  x„  =  x„o  +  6xn,  the  difference  equation  for  the  error,  8: r,  in 
the  parallel  direction  is 


6(i!  -  2x0  +  £-i) 


(6xo  •  V)-9^A<2 


m 


(±V  ■  EAt2)6x0 
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—4v6ne2At2 

—  - ox  o 

m 

^  _  _ 

The  maximum  8E  occurs  for  =  N8x0  or  8n  =  n0.  This  is  equiv¬ 

alent  to  displacing  all  the  particles  together.  We  may  obtain  a  dispersion 
relation  for  this  scheme  by  (z)  transforming  to  8xn  =  8xQe%nuAt  [21,  page 
664], 

eiu,At  _  2  +  e-i^t  =  _(WpeA<)2 
—4  sin2  u  At/2  = 

requiring  u>peAt  <  2  for  stability  to  parallel  displacements.  A  more  careful 
analysis  reveals  that  the  actual  stability  requirement  is  wpeA t  <  1.62  [5,  p. 
184]. 


A. 2  Gyrokinetic  Solution 

A  similar  analysis  shows  perpendicular  displacements  to  be  stable.  How¬ 
ever,  the  full  dispersion  relation  includes  non-physical  effects  from  aliasing 
in  frequency  which  causes  instabilities  analogous  to  cyclotron  waves  [5,  p. 
201].  In  general,  aliasing  can  be  reduced  by  choosing  smaller  values  of  ft  A  t 
[5,  p.202-203].  Thus,  in  particle  pushing  schemes  including  gyromotion  the 
time  step  is  limited  by  the  time  to  complete  a  gyrorbit.  It  is  desirable  to 
use  a  scheme  which  eliminates  this  stability  restriction. 

Such  a  scheme  results  from  an  (analytic)  averaging  over  particle  orbits. 
When  we  average  the  single  particle  equation  of  motion  over  the  gyromo¬ 
tion,  to  order  e  ~  —  ~  ^  ~  p/L  (with  L  =  \B/VB\  [31]), 

R  =  -[E(R )  +  (k/c)  x  B(R)]  -  (^-)VB(R)  +  o(e) 
m  m 

[32,  eqn  (12), p.  84]  where  p  =  This  equation  is  as  difficult  to  solve 
as  the  original  equation  and  carries  unphysical  properties  with  it  [32].  A 
substitute  for  this  equation  is  the  coupled  set,  to  o(e2), 


and 


f?  e“i  .  e  -*  p  ~ .  dej  dun,. 


[32,  eqn  (17),(20)]In  our  model  B  is  stationary  and  both  £  and  i?  are  in 
the  p,  z  plane  so  that  U£  =  U£<£.  Invoking  azimuthal  symmetry  =  0) 
the  parallel  equation  reduces  to 


duji 

dt 


eEw 


m 


n  dB  -  - 

+  uE  •  ( uE  •  V)6 


m 


5s 


and 


k  _  r  ^X  p(V£)X  vj|,56  *>||  due  bz  2i7 

^x  =  [-c~g-+  — . +  7^(^:)x  +  7T(— K  “  77 Zue]4> 


mft 


ftv5s 


ft  5s 


+  ( 


ftp 

1  5lt£ 


ft  5f 


+ 


UEbpV\\ 

pQ, 


)(b  x  *) 


where  X  denotes  the  direction  (j>  x  5.  All  except  the  last  two  terms  in  R± 
are  in  the  <f>  direction.  In  general,  these  equations  are  not  amenable  to 
time-centering. 

Formally  these  equations  were  derived  for  ~  e,  E_l  ~  1.  In  some 
instances  the  additional  assumption  that  E±  ~  E\\  holds  and  the  last  terms 
smd  others  involving  ue  may  be  ignored.  In  our  simulation  we  assume 
instead  ^  ~  u>  ft,  [o(^)  ~  e2],  where  u>  is  the  time  scale  for  large  field 
changes  so  that  the  first  of  these  terms  may  be  ignored.  Additionally,  for 
weak  curvature  the  second  term  is  on  the  order  of  p/r  ~  e  in  comparison 
to  other  terms  in  the  perpendicular  velocity. 

Writing  the  velocity  in  cylindrical  coordinates,  we  have 


v  =  rr  +  r<j><j>  +  zz 


and  the  parallel  and  perpendicular  velocities  separate  into  their  components 


V|, 


V±  = 


rr  +  zz 


r<f><j>  =  Vd  + 


q  (E)  X  b 
m  ft 


where  vd  =  +  ?f(f!)x  +  Tf(^)x  -  fcu2E]j>  is  the  <t>  component  of 

the  perpendicular  velocity  less  the  E  x  B  drift. 
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Differentiating  the  parallel  velocity 


dv< ||  _  d(u„6) 
dt  dt 
db 
dt 


rr  +  zz 


r i  -  r_i 


At 


—r  + 


z  i  —  z_i 


At 


dv  |i  j  d6  „  .x  .x 

—7^-6  +  vii—  =  rr  +  zz  +  rr  +  zz 
dt  "  dt 

d  .... 

— (6rr  +  6*2)  =  bTr  -f  6rf  +  6*2  +  6*2 
at 

^6  +  i>||(6rr  +  6*2 ) 

dv\\l,  2,dbr*^dbz 

-dTb  +  v^-dIr  +  ^ 

dv  in 


ds 


(1  +  cos2  9)bg  « 


<A3> 


where  we  have  used  t?n  •  z  =  U||,6*  =  2.  We  observe  that  the  above  equation 
consists  of  an  acceleration  and  a  pure  rotation  so  that  we  may  apply  the 
Boris  algorithm  with  v±  — ►  t>||  and  ft  — ►  3ujj(l  +  cos2  9)bg/r(l  +  3 cos2  9). 
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Appendix  B 
Poisson  Solvers 


B.l  Charge  Aggregation  and  Force  Interpo¬ 
lation 

B.l.l  Pairwise  Force  Calculation 

The  most  obvious  way  to  calculate  the  force  on  a  particle  is  to  calculate 
the  force  from  the  local  B  field  and  the  resultant  pairwise  electric  field. 
Naively  this  method  involves  o(n2)  ~  calculations  for  n  particles. 

In  one  dimension  the  number  of  calculations  may  be  reduced  to  ~  n  log  n  by 
particle  sorting.  In  higher  dimensions  we  are  unable  to  use  this  technique, 
and  even  in  ID  the  number  of  calculations  may  be  prohibitive. 

B.1.2  E  Field  Grid 

At  the  expense  of  some  loss  in  precision  an  alternate  route  may  be  taken. 
A  grid  may  be  superimposed  upon  the  model,  then  charges  are  assigned  to 
and  Poisson’s  equation  is  solved  at  each  grid  point  [23].  Besides  avoiding 
computation  of  pairwise  forces,  the  grid  facilitates  study  of  field  behaviour 
in  space  and  time  and  smooths  collisions  among  the  limited  number  of 
particles. 
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B.1.3  Direct  Integration 

Once  we  know  the  electric  field  at  grid  points,  we  must  calculate  the  force 
on  particles  between  points.  We  could  interpolate,  using  a  Taylor  series 
expansion.  An  equivalent  approach  in  one  dimension,  with  strong  physical 
appeal,  is  direct  integration. 

f dx(V  ■  E)  =  E(x)  -Eg  =  4ir  [* pdx 

JXg  JXg 

[x9+ 1 

Similarly  Eg+ %  —  Eg  =  pdx  =  47r  <  p>  A 

Jig 

or  E(x)  =  Eg  +  (Eg+1  -  Eg)6  -  4tt  [  (<p>-p)dx 

JXg 

=  6Eg+l+(l-6)Eg  +  o(A 2) 

where  6  =  and  <  •  •  •  >  denotes  the  average.  Thus  the  particle  force 

is 

F(x)  =  ,[<£,+1  +  (1  -  «)E,] 

While  it  is  possible  to  mix  schemes  for  charge  sharing  and  force  cal¬ 
culation,  due  to  the  reciprocal  nature  of  fields  and  density,  use  of  a  single 
scheme  for  both  prevents  unphysical  self  forces  from  arising  [5,  p.162-3]. 


B.2  Electric  Field  Calculation 

B.2.1  Difference  Equations 

In  the  ID  case  we  may  integrate  Poisson’s  equation  (3.1(a)),  suitably  in¬ 
terpolating  p  between  grid  points.  This  approach  is  equivalent  to  a  finite 
difference  solution, 


E9  ~  E9+ 1  =  Pg+  1&X  +  o(Al2), 

where  we  have  used  the  difference  analog  to  the  operator  Ag+^  —  Ag_  i, 
A  a  quantity  defined  at  the  grid  point  g. 

In  the  ID  electrostatic  case  Poisson’s  Equation  (3.1(a))  becomes 

-  2 <t>g  +  4>g- 1  =  -4irpgA2 
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We  may  obtain  the  solution  for  <j>g  after  inverting  the  resulting  tridiago¬ 
nal  equation.  The  Buneman  algorithm,  as  presented  in  Appendix  C,  is  a 
particularly  efficient  application  of  this  scheme  in  more  than  one  dimension. 


B.2.2  Fourier  Transform  Solution 

Discrete  Fourier  Transfer'  is 


Solution  of  Poisson’s  equation  can  be  obtained  by  using  discrete  Fourier 
transforms,  just  as  one  would  use  Fourier  transforms  in  the  continuous 
sense.  The  discrete  Fourier  transform  of  (B.4)  on  the  interval  x  =  0,  L, 
with  grid  points  at  xg  =  <7 A  +  xo,  g  =  0, 1, 2,  •  •  • ,  N,  is  obtained  from  the 
continuous  (finite)  transform  as 

Pk  =  jL  dx'p(x')e-ikx' 

Jo 

=  Eh  ja  dx'6(x'-x ,)e~,kx' 

=  *Ewik“ 

3 

The  inverse  is  found  by  noting  that 

eik(xg-xg,)  __  (  N  9  =  9  +  IN 
\  0  g^g'  +  lN, 


% 

E 


l  an  integer,  so  that  pg  =  ^  Ylk  Pkeikx 
m  =  l,2,---,JV  and  Afc  =  f. 


*,  where  k  takes  on  the  values  mAh, 


To  solve  a  linear  equation  such  as  Poisson’s,  we  apply  the  operator 
AJ2ge~'kx',  obtaining  the  k-space  solution  as 


£9  =  f 
gk  =  Cllfk 

We  then  apply  the  inverse  transform  to  obtain  g.  On  its  face  this  procedure 
appears  to  require  ~  N2  operations  (~  N  summations  for  each  k  x  N  ks.) 
However,  the  Fast  Fourier  Transform  (FFT),  has  decreased  that  number  to 
~  N  log  iV,  a  huge  savings[16].  We  shall  write  the  FFT  algorithm  as 

FFT(A)  =  n,m  discrete 
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We  are  able  to  compute  the  transform  of  quantities  defined  at  points  dis¬ 
placed  from  the  grid  by  including  a  phase  factor: 

£  e~'kx°Ag  =  £  e~ik{9A+s)Ag  x=gA  +  6 
9  a 

= 

9 

=  FFT(Age~ikS) 

Ak  =  e~ikSFFT(Ag) 

Similarly,  the  phase  factor  must  be  included  when  we  find  the  inverse  trans¬ 
form 

These  notions  are  useful  when  examing  charge  sharing  and  interpolation 
schemes,  but  for  grid  quantities  where  x=xg>  =  g'A  +  6,  we  have 

A(xg)  =  53  e22^21  FFT(Ag)  =  FFT~1[FFT(Ag)\ 

k 

and  the  phase  factor  is  unneeded. 


ID  Green’s  Function  Solution  to  Poisson’s  Equation  with  Mixed 
boundary  conditions 

Continuous  Solution  We  desire  to  solve  the  continuous  Poisson’s  equa¬ 
tion  for  E  and  use  it  as  a  basis  for  the  discrete  solution  in  our  simulation. 
Thus,  we  seek  the  solution  to  the  ID  boundary  value  problem: 


-V2<£  =  47rp,  «K0)  =  0,  <£'(/)=() 


The  solution  may  be  expressed  in  terms  of  a  Green’s  function  which  satisfies 
the  equation 

V^G  =  -4tt£(x  -  x')/A,  (B.l) 

with  boundary  conditions  on  G  to  be  determined.  Defining  g  =  fA  Gda\ 
we  multiply  (B.l)  by  <£  and  integrate  by  parts  to  obtain 


<£(x)  =  p(x')G(x,x')dV  +  ~  ^[G(x,*')~  -  <f>{x')  — — ~ ] da1 
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n  =  -x 


A  A 

n*  x 


Figure  B.l:  This  sketch  shows  the  area  of  integration  used  in  evaluating 
the  Green’s  function  solution  to  Poisson’s  equation. 

=  jf  <M*)»  +  -  ^^)lot[24,p.41] 

At  x=0  where  0  =  0  it  is  most  convenient  to  remove  the  second  term  by 
requiring  G=0.  At  x=L  where  <f>'  =  0  we  can,  at  best,  satisfy  the  boundary 
condition  imposed  by  Poisson’s  Equation: 


Jf&G-fgut  -  -4, 

d9  |  dg  _ 

-  ~4jt 


This  implies 


<f>(x)  =  Jv  p(x)G(x,  x')<Px'  +  <f>(L)[  1  - 


Evaluating  this  expression  at  x  =  0  with  the  boundary  condition  G(0,  x')  = 
0  gives  |f  |0  =  4tt  or  |f  | l  =  0  just  as  for  the  potential  <f>. 

We  solve  for  G  using  an  alternative  Green’s  function, 


E( x)  =  f  dx'd(x  —  x')p(x'), 
Jo 


where  d(x  —  x')  is  the  solution  to 


and  g  and  d  axe  related  by 


-(d)  =  47r5(x  —  x') 
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Applying  the  condition,  E(L )  =  0, 


d{L  -  x')  -  d(x  -  x')  =  [L  4irS(x"  -  x')dx" 

Jx 

1/  /\  f  0  x  x 

v  '  [  4ir  i  <  i 


After  integrating  and  applying  the  condition,  <^>(0)  =  0, 


g(x  —  x')  —  g{ 0  —  x')  =  —  [  dx"d(x"  —  x ') 

Jo 

{4ttx'  x  >  x' 

4irx  x  <  x' 


Discrete  Solution  Now  we  may  solve  for  E.  Making  use  of  a  continuous 
finite  transform  in  x  and  a  discrete  transform  in  k,  we  have 

£(x)  =  [Ldx\  !>«•'*'*') 

Lk  w  Jo 

=  7E<W‘* 

L  k 

or  Ek  =  dkpk 

In  other  words,  to  solve  for  the  electric  field  of  some  arbitrary  charge  dis¬ 
tribution,  we  need  only  solve  once  for  dk  and  algebraically  find  the  solution 
for  the  particular  charge  spectrum. 

Proceeding  to  the  solution  for  dk ,  we  begin  by  solving  for  the  discrete 

*• 

[L 

<f>(x)  =  I  dx'g(x  —  x')p(x') 

Jo 

=  f  dx’g(x  -  x')(j'52pke,kx') 

JO  L  k 

=  VE  M  /  dx'x'  +  /  dx'x)etkx>  +  p0(  f  dx'x'  +  f  dx'x)\ 

±J  LJrt  JO  Jx  Jo  Jx 
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E  =  —V<t>  =  ^E-^r^-^  +  Mx-i)]  (B.2) 

These,  of  course,  may  be  recognized  as  a  particular  solution  to  the  inhomo¬ 
geneous  equation  (the  periodic  (in  x)  parts)  and  a  general  solution  to  the 
homogeneous  equations  (the  parts  non-periodic  in  x.) 

An  “Exact”  Poisson  Solver 

In  many  instances  use  of  the  Fourier  Transform  solution  is  preferable  to 
that  of  difference  equations  which  results  from  approximating  the  operator 
gj  (or  -£s)  when  solving  for  <f>.  As  its  name  implies,  the  FFT  is  relatively 
fast,  particularly  in  multiple  dimensions  and  unlike  other  fast  solvers  such 
as  Buneman’s,  allows  easy  analysis  of  the  field  spectra,  and  incorporation  of 
particle  shaping.  If  we  plan  to  use  FFTs  to  solve  these  equations,  however, 
we  must  be  prepared  to  live  with  aliasing  which  result  from  inability  to 
distinguish  short  wavelength  modes  from  others  on  a  finite  grid.  Indeed, 
we  can  readily  show  examples  using  Fourier  Transforms  where,  despite  the 
advertised  (local)  error  estimates,  global  errors  can  be  disastrous. 

In  the  previous  section  we  implicitly  solved  for  d*  =  —  ^  in  equations 
B.2  and  B.2.  In  general,  however,  the  choice  of  the  Poisson  solver  is  some¬ 
what  arbitrary.  Eastwood  [19]  makes  a  case  for  choosing  the  (periodic) 
solution  to  Poisson’s  equation  as 

Ek  =  dkpk  {B.  2') 
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where  c4  minimizes  the  mean-squared  deviation  of  F,  the  resulting  force, 
from  R,  a  “reference”  force.  In  our  case  the  reference  force  is  that  defined 
by  the  particle  shape  and  E  field  (B.2). 

The  optimal  <4  =  |£4|2  =  Rk/\  |£4„|2|2.  Here  Rk  is  the  Fourier 
transform  of  R  (with  the  assumption  that  R  is  band  limited),  £4  is  the 
transform  of  the  shape  function  U,  and  £4n  is  the  alias  of  £4  due  to  finite 
transforms.  Basically  <4  is  that  from  the  Green’s  function  solution  to  the 
Poisson’s  equation  adjusted  to  account  for  particle  shape/force  weighting 
schemes  and  aliasing  due  to  finite  grids. 

The  assumption  of  a  band-limited  force  is  tenuous  from  the  outset,  and 
the  proper  choice  of  Ek  for  non-periodic  boundary  conditions  is  difficult  if 
not  completely  ambiguous.  Thus,  in  many  cases  the  direct  k-space  solution 
to  Poisson’s  equation  (B.2)  is  likely  to  be  as  good  as  (or  better)  than  other 
more  sophisticated  schemes  [19,  p.15]. 

Trying  to  retain  the  advantages  of  the  FFT  without  loss  of  accuracy  due 
to  aliasing,  we  take  a  slightly  different  approach.  Rather  than  approximate 
the  operator  and  take  the  charge  distribution  as  known  exactly  at  points,  we 
shall  use  the  exact  operator  and  approximate  the  charge  distribution  over 
the  continuous  interval  about  the  known  points,  xg  ~  gA  +  6,  g  =  0, 1, 2  •  •  •. 
We  take  as  given,  pg,  and  expand  the  continuous  function  p  about  these 
values.  Of  course,  if  we  know  p(x)  everywhere,  we  could  (conceivably)  solve 
for  E  exactly. 

In  one  dimension 

dE  .  .  .  rrV-,  dp.  d2p.  (x  —  Xg)2  .  . 

qx  ~  4?r p(x)  —  47r{[£]  p(xa)+  —  |x=Ij(x— xg)+-^^\x-Xg  —  I 


where  xg  is  a  grid  point,  A  the  grid  spacing,  and 


n 


9 


1  |z-xa|<! 

0  jx-x3|>f 


is  the  “hat”  function.  Further  we  may  require  that,  if  we  find  the  value  of  pg 

by  “adding  up”  particle  charges  within  each  gth  cell,  fg&+s'll  pdx  =  pa- 

This  effectively  eliminates  all  even  values  of  the  expansion  (except  for  the 
zeroth  term.) 

Now,  to  solve  Poisson’s  equation,  we  use  Fourier  transforms  .  However, 
we  would  like  to  avoid  the  problems  of  aliasing  or,  alternatively,  problems 
of  truncation  associated  with  difference  schemes. 
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Specifically,  the  FFT  and  continuous  transforms  for  p  can  be  shown  to 
be  related.  Applying  the  operator  /0L  dxe~tkx  •  •  •,  we  obtain  the  transform 
of  the  quantities. 

_  [L  -t'fcx  r— i  rV''  -  I  I  /_  _  \  I  1 


=  Jo  dxe-ik*  ng  E  pg  +  -£\a{x  -  xg)  +  . .  •] 

rgA+S+y  Qn 

=  dx^+Q^\3^-^9)  +  --^  'X 


We  may  evaluate  the  integral  piece  by  piece  for  each  cell,  xg  -  ~  <  x  < 


xg  +  and,  taking  pg,  ||  |a  •  •  •  as  known  grid  quantities,  we  find 
r9+i  dxe~ikx  =  )  _  e-«*(*9-f-)l 

Jxn-%  ~IK 


t  .‘k><  /  v 

=  — e  9(e  ’2  —  e  *  ) 

K 

*  ttx  ■  ,kA..,kA. 
=  Ae  ,te*sin(— )/(—), 


a  well-known  result  [5,  p.  169]. 


J*3\7  dxe~ikx(x  -  xa)  =  J3^  dx'e-ik(x'+x^x' 

=  e~ik*9  / 1  x'dx'e-ikx' 


-ikxa:_d  f  1  -ilex' 


=  e  — n 


-  P  dx'e 

dkJ 


-ilex  ■  d  r  •  .fcA.-.fcA .  1  . 

=  e  ”5*lsm(— : '/(T)1a 


=  c-,fcx»i[— [cos(-— )  — 


T^vc-f) 

T  z 


sin(¥) 


=  e  ’**9-[cos(— )- 


Assuming  5  =  0, 
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=  £[j^"‘/2  -  ‘“'V'*"  +  J_\  dAxe^’] 

...  Y--rsin(¥)  ,  .  9p..d  sin(fcA/2)  „ 

=  £[/>„/»(*)  +  ^U.(t)  +  •  ■  ■]e-ta>A 


where/„(fc)  = 

Note  that,  although  we  do  not  have  a  measure  of  dp/dx\g  and  higher 
derivatives,  we  may  either 

1)  approximate  their  values  with  a  difference  equation,  for  ex¬ 
ample,  or 

2)  use  their  exact  values  in  terms  of  the  unknown  /?*. 

We  take  this  second  approach  for  now  but  return  to  the  first. 

Noting  that  p  =  jr  YlkL-<x>  e'kx  Pk,  we  have 

d  1  00 

(5jm  =  T  £  •  •  • 

k=—oo 


Performing  the  operation  A  J2g  =  LSkk1,  f>kk>  the  Kronecker  delta, 

x  _  /  1  k  =  k' 

6kk'  ~  {  0  k^  k ', 

we  have 

Pk  =  £A,/o(*)e-ik*'  +  £[^(*'4)”/oWW 

9  nodd  n •  aK 

We  recognize  the  sum  over  n  as  (formally)  [sin(-k  ^)f0(k)]  where  sin(-fc^) 
operates  on  fo(k).  We  finally  have 

„  -  ,  f  /<>(*)  1 

k  m  1  -  sm(-k^)fo(k) 

That  is,  the  infinite  transform  pk  is  related  to  the  finite  transform  pm  by 
the  relation 

Pk  =  PmF(k) 
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This  is  a  pleasant  state  indeed.  Let  us  assume  that  we  desire  to  obtain 
the  function  Eg( say)  =  J2V=-<x>  e>kx3Ek-  We  know  that  Ek  and  Em  are 
related  through  the  aliases  to  Em,  that  is 

oo  N/2  oo 

E,=  £  e‘k’<Et  =  £  £  e<k>x>Eh 

k=- oo  m=~N/2+l  p=-oo 

where  kp  =  m  +  Np.  If  we  can  solve  analytically  for  £*(e.g.  Ek  =  —  ^r4-), 
then  we  may  solve  “exactly” 

OO 

E,  =  Y, 

m  p=— oo 

=  E  W/M 

m  p=-oo 

m  m  p  1 -sin(-fcp^)/o(A:p)J 

__  «mxa/  -Pm  X  y*(  -JT  \ _ (  — l)P/o(m) _ 

We  recognize  the  terms  in  the  solution: 

1)  Em  e,mi»(-i^)  is  the  standard  FFT  solution  for  E  at  the 
grid  points  without  any  corrections  for  aliasing. 

2)  5Zp  •  •  •  is  the  correction  factor  for  aliasing. 

3)  sm(—kp-^)fo(kp)  carries  the  degree  to  which  the  continuous 
p  is  approximated  by  the  Taylor  series  expansion. 

By  truncating  this  series,  we  are  (equivalently)  creating  a  difference  equa¬ 
tion  to  solve  for  Eg.  This  procedure  is  equivalent  to  creating  a  band- limited 
function  as  posited  by  Eastwood. 

We  illustrate  this  point  by  dealing  with  only  the  zeroth  order  term  in 
the  expansion  for  p{ i.e.  p  =  l~l gpg ). 


E,  =  JV-.  (-,•£)£ 
m 


x  +  up 
m  ml 


(-l)P/o(m) 
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where  we  have  used  the  expansion  [10,  p.  217] 


cot  x  = 


E 


i 

X  +  TTTl 


We  are  not  surprised  that,  in  agreement  with  Eastwood,  this  solution  cor¬ 
responds  to 

<f>g+ 1  ~  2<Ag  +  <f>g-\ 

P9  ~  A2 

„  _  <f>g- fl  ~  4>g- 1  , 

Eg  2A 

We  may  obtain  as  accurate  an  expression  as  patience  allows  by  evaluat¬ 
ing  the  sine  series  to  any  desired  order  and  physically  summing  the  resulting 
expression  over  as  many  ps  as  we  choose!  If  we  intend  to  use  this  scheme 
in  a  program  where  we  repeatedly  evaluate  Ea  given  pg,  we  need  only  eval¬ 
uate  the  modifying  terms  at  the  start  of  the  program  once.  Since  this  term 
converges  for  p  as  we  expect  rapid  improvement  in  our  evaluation  for 
E  over  its  aliases.  We  may  in  turn  increase  our  accuracy  in  terms  of  A” 
to  any  compatible  order  by  keeping  higher  and  higher  terms  in  the  sine 
expansion.  For  example,  keeping  the  first  sine  term  should  give  an  order  of 
accuracy  better  results  than  the  common  finite  difference  scheme. 

Furthermore  we  still  retain  the  ability  to  use  particle  shaping  in  these 
schemes.  The  subject  of  particle  shaping  remains  a  separate  issue  -  what 
particle  shape  gives  the  best  physics  shall  be  briefly  considered  in  the  next 
section. 


B.3  Particle  Shapes 

B.3.1  Cloud-in-Cell  Weighting 

At  each  time  step  the  particle  contributes  to  the  density  and  derives  a  force 
from  the  fields  based  on  its  position.  Commonly  the  rule  for  assigning  the 
density  and  calculating  the  forces  for  a  discrete  array  of  grid  points  based 
on  a  continuous  particle  position  is  known  as  interpolation.  In  particle 
simulation,  however,  it  is  useful  to  consider  the  alternate  view  of  finite  size 
particles. 
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In  our  ID  case,  to  calculate  the  force,  we  assign  particle  charge  to  each 
of  the  grid  points: 


Qg  = 


! 


?(!  -  kiEii) 
0 


|x,'  Xg\  <  A 

l*i  -  Xa\  >  A 


(B-3) 


The  same  scheme  is  obtained  by  considering  the  particle  to  be  a  constant 
density  slab  of  thickness  A  and  infinite  cross  section.  This  is  known  as 
cloud-in-cell  (CIC)  weighting  [5]. 

Charge  density  for  field  calculations  may  also  be  assigned  to  each  grid 
point  based  on  the  volume  of  this  particle  in  each  cell.  More  generally,  after 
counting  particles  within  each  grid  cell,  we  obtain  a  charge  density, 

Jj  A  A 

p(x)  =  Z  /  dx'p(x')  j  ’  dAxS(x'  -xg-  Ax)/  j*  dAx  =  £  qa&(x  ~  xg) 

g  Jo  2  J  2  g 

(BA) 

where  S  is  the  particle  shape  which  determines  assignment  of  particle  charge 
to  the  grid  for  calculation  of  field  quantities. 

The  idea  of  finite  size  particles  is  easily  incorporated  into  plasma  theory 
and  is  particularly  useful  in  studying  the  effects  of  higher  order  interpola¬ 
tion  schemes  or  the  use  of  “smoothing  factors”  to  promote  more  realistic 
simulation  behaviour.  For  example,  the  use  of  “smoothing  factors”  may  be 
interpreted  in  terms  of  “effective  particle  shapes”.  Thus,  a  principal  benefit 
of  the  finite  size  particle  concept  is  as  a  bookkeeping  tool.  It  ensures  that 
we  may  analytically  compute  the  effects  of  our  charge  sharing  and  force 
interpolation  schemes. 


B.3.2  Optimal  Particle  Shape 

Potential  for  Finite  Size  Particles 


The  Coulomb  potential  for  a  pair  of  point  particles  is 

<t>(r)  =  — ,  r  =  |x2  -  fi| 
r 

This  can  be  written  in  the  equivalent  form 


<f>  =  J  dx'  J 


dx 


„e2S(x'  -  xj )6(x"  -  x2) 
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for  point  particles  or,  more  generally, 

jl  f  f  ~  X\)S{x"  -  X2) 


4>  =  J  dx'  J  i 


for  finite  size  particles  which  includes  point  particles  as  a  subclass. 

To  exhibit  the  explicit  dependence  of  p  on  the  separation  of  finite  size 
particle  centers  we  note 

S(x' -  x,)  =  (±)3  J dkSkeikl*'-*'' 

and  therefore 

e2  r  j  >  r  ,  » r r  ,l,sks'ke,k^x'-x')cik',r"-x^ 

(2tt)6  J  J  J  J  \x"-x'\ 

Changing  coordinates  to  x  =  x  |y  ,  the  position  of  the  center  of  charge 
elements,  and  Ax,  the  separation  between  charge  elements,  we  have 


e2  f  j _  f  j  a  f  j,  f  ,SkSk>tik(?-*f"x^ik'(x+*r-x^ 

*  =  iwJJ  JJ  - W\ - 

=  (£rhldA*Jdk  I 

We  next  use  the  identity  /  dxe~^k~k'^x  =  f  dxe~,kx(e,k"x)  =  2i rS(k  —  k'). 
(This  can  be  most  easily  seen  in  that  (j^)J  dke,kx6(k  —  k' )  =  e,k'x' /2n  so 
that  the  Fourier  Transform  of  e,k'x  is  2ir6(k  —  k ')). 


6  -  6 2 
^  (2tt)6  j 

1  dAx  j 

(  dk 

<0 

N  fc. 

II 

1  dAx  j 

1  dk 

Jdk" 

e2 

(2tt)3J 

f  dAx  j 

f  dk 

SkS-k 

|Ax| 

Now  we  may  transform  <f>{r) 

4>k>  =  J  dre-ik'r<j>(r) 
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Two  Body  Corellation  Function  for  Finite  Size  Particles 

The  use  of  finite  size  particles  to  reduce  collisions  in  simulations  is  a  common 
technique  to  compensate  for  reduced  shielding  of  charges  from  the  relative 
smallness  of  the  plasma  parameter  -  nAf,.  We  wish  to  examine  how  the 
two  body  correlation  function  (and  therefore  collisions)  behaves  with  finite 
size  particles  and  postulate  the  existence  of  an  “optimal”  particle  shape  in 
regard  to  the  correlation  function. 

In  thermal  equilibrium  Frieman  and  Book  [20]  have  shown  that,  assum¬ 
ing  a  correlation  function  of  the  form 

-  zal'fi't*)  =  /m(vi)/m(u2)^(|3:i  -  i2|) 

where  }m  =  (m/2TrkT)$e~m'^f2kT,  the  spatial  dependent  part  of  ^  obeys 
the  equation 

.  .  d*  1  d<f>  ,  .  T  l  d<j>  . 

or  ("■  - "»>  •  iff + +  ^  =  J  -  **\m&*\) 
jp  +  wi£(*  +  1)  = 

The  effect  of  finite  particle  size  on  the  potential  <j>  is  to  reduce  it  at  short 
distances  so  that  it  remains  finite  even  at  zero  separation.  We  plot  the  elec¬ 
tric  field  and  potential  versus  separation  distance  in  the  figure  below.  We 
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Figure  B.2:  This  shows  the  potential  vs  radius  for  particles. 


note  that,  as  the  particles  overlap,  the  force  decreases,  finally  approaching 
zero  for  small  r.  At  large  distances  the  electric  field  and  potential  both 
must  approach  the  classical  Coulomb  potential. 

If  is  on  the  order  of  e  at  a  (by  assumption)  then  it  is  of  the  same 
order  for  r  <  a.  This  then  changes  our  assymptotic  analysis.  The  terms 
above  are  in  the  ratios  given  below,  as  deduced  from  Book  and  Frieman, 

f*  f  nr03fc/rfA*V(|*-A*|)¥(|Ax|) 
“jinairro  <  n~i  e  e  <  e2 

r0  >  a  ~  Ad  e  e  e 


where  V  =  x  =  — . 

K1  7  VQ 

For  “small”  r  we  have  the  balance 


d*  dVfr 

dx  +  ai^  +  1) 

vd(*  +  l)  dev,T 


dx 


=  0 


=  0 


[ev(*  +  l)]  =  0 


¥ 


— 1  +  conste  v 
— 1  +  conste 


Note  for  V  large,  balances  |j,  while  for  V  small,  effectively  balances 
itself. 
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For  “large”  r  (V  small)  the  balance  must  be 
flvlr  f)V  dr 

fc  +  fc  +»*iS/^V(|£-  AxlWIAxI)  =  0 

We  have  immediately  that 

'£  +  V  +  nA^,  J  dAxV(\x  —  Axl)'I'(IAxl)  =  const 
Using  the  convolution  theorem, 

J dAxV(\x-  Ax|)*(|Ax|)  =  J dkVkVkeikx, 
we  have  the  solution  in  Fourier  components 

Vk  +  Vk  +  nVk'&k  =  [47 r2  8(k)8(n)  l  k2)consi 


1  +  nl4  A:2  1  +  n!4 


*  =  _(J_)3  /•<,*[ - ^ - ,to 

V2tt'  J  ll  +  nV*  A:2  1  +  nV*  J 

=  -  _L  [  dkk2r  Vk-[(eikr  ~  e— ,fcr)/*A:r] 

1  +  nVo  47T2  J  1  +  nl4 ’  J 

For  <f>  we  have  the  standard  Coulomb  potential  integrated  over  the  shape 
functions 

V  =  /  dAxx  j  dAx2— e^Axi)5(Ax^- . /ArT 

./  ./  |xi  +  Axi  —  x2  —  Ax2| 

th  =  j  dxje^ 

=  27 Te*  J  drr  J  dfie~,kr,t 

=  ^  r  dr{e~ikT  ~eikr) 

—ik  Jo 

2ne 2  A-ik-a)r  (i \k-a)r 

=  ~r~i  lim[— - - ]IS° 

A:  cr-.o  — —  a  tA:  —  a  J  0 
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so  that  *  =  /“  /(fcAD)a](^^). 

For  Si  =  1  (point  particles)  we  have  poles  at  k  —  ±i/ \p.  We  extend 
the  limits  of  integration  from  —  oo  to  oo  since  the  integrand  is  even  in  k  and 

1  e2  -oo  U2  ikr  _  -ikr 

*  =  kT  I- oo  ^Jb2-f  (1/AD)2^  ifo* 


We  note  that  the  portion  of  the  integral  involving  e,kr  must  be  evaluated 
over  contour  1  while  that  for  e~,kr  must  be  evaluated  over  contour  2  in 
Figure  B.3. 


dk-^—C- 

k 2  +  (1/Ad)2  ikr 


iff 


<(*)-• 


3 

1  /2ir 

r  dei - 


Jo 


o-r/^D 


With  a  like  contribution  from  the  other  contour  so  that 


*  = 


kTr 


rA  d 


as  in  Frieman  and  Book[20]. 


An  “Optimal”  Particle  Shape 

We  see  that  the  trade-off  between  the  (real)  point  particle  case  and  the 
finite  particle  case  must  be  the  density  (in  the  form  of  Ad  =  yT/n)  and 
the  particle  shape.  Since  g  (the  two  body  correlation  function)  represents 
the  error  due  to  collisions,  we  should  seek  to  minimize  this  error. 

An  appropriate  shape  should  result  from  the  minimization  principle: 

A/rf>x(*w -»(«))' =  o 
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B.4  Extension  to  Two  Dimensions 


We  extend  the  simulation  to  two  dimensions.  The  model  is  similar  to  the  ID 
case  but  now  the  potential  varies  both  in  r  and  z.  The  model  and  boundary 
conditions  are  pictured  in  Figure  6.1.  We  shall  find  the  boundary  condition 
at  r=R  is  flexible.  The  condition  at  r=0  reflects  the  fact  that  the  particle 
distribution  is  continuous  (no  line  charge). 

B.4.1  Poisson’s  Equation  in  Spherical  and  Cylindrical 
Coordinates 

Azimuthal  symmetry  about  the  z-axis  suggests  use  of  cylindrical  or  spheri¬ 
cal  coordinates.  Poisson’s  equation  and  the  electrostatic  potential  in  these 
two  systems  is: 

or 

1  3  .  2  _  .  1 

-TTr^r  +  — ——(sm6E6)  =  4t rp 
r2  or  r  sin  0  30 

E  dr  +ra/ 

Numerical  solution  of  Poisson’s  equation  is  facilitated  when  the  V2  operator 
can  be  expressed  with  one  “Poisson-like”  variable.  Cylindrical  coordinates 
therefore  are  chosen. 


B.4. 2  Continuous  Solution 

In  cylindrical  coordinates,  we  have  an  equation  for  the  potential  <f>(p,  z): 

ld(dt  d*<f> 

-pTp(pp  +d?  =  -4Kp 

The  charge  density  p  can  be  expanded  in  an  infinite  Fourier  series,  £  YlT=-oo  p(p)  exp,fe% 
with  pk{p)  =  fo  dzexp~'k*  p(p,  z ).  This  suggests  a  form  for  the  (f>  solution, 

^  =  7  £  4>k{p)explkz +<i>h, 

L  k=- oo 
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where  </>k(p)  are  particular  solutions  to  the  inhomogenous  equations 

1  d  d(f>k{p )  iiit  \  i  t  ^ 

pTpil’-ar)  - k  Mp)  =  -4kpM 

and  4>h  is  a  homogeneous  solution  to  the  cylindrical  Poisson  equation  which 
allows  <f>  to  satisfy  the  boundary  conditions  in  z,  <^(0)  =  0 ,<f>'(L)  —  0. 

The  Green’s  function  solution  for  <fik  may  be  expressed  in  terms  of  mod¬ 
ified  Bessel  functions  Io(p),  A'o(p).  The  requirement  for  finite  <f>  at  r=0 
precludes  use  of  Kv  which  behaves  as  ln(|)  for  i<Cl. 

We  may  add  to  this  any  solution  to  the  homogeneous  equation  which 
satisfies  the  given  boundary  condition.  However,  we  choose  this  solution 
to  be  zero  in  absence  of  a  source,  p.  The  homogeneous  solution  consists 
of  products  of  exponentials  (real  or  imaginary)  in  z  and  appropriate  Bessel 
functions,  <f>h  =  7 Zk>  Jo{k'  p)[Ak-ekz  +  Bk<e~k  *].  The  boundary  conditions  on 
the  homogeneous  solution  derive  from  those  on  <j>  itself,  i.e.  <f>h  =  <j>  —  <f>p. 
The  boundary  conditions  in  r  are  $,(0,  z)  =  0  and  <j>h{L,  z)  =  0.  The  first 
condition  allows  solutions  of  the  form  Jo  while  the  boundary  conditions 
at  r=R  restrict  kR  to  zeroes  of  J0.  The  boundary  conditions  in  z  axe 
MO)  =  ~T,k  4>k{p)  and  ML)  =  ~i  72k  Hk{p)  and,  using  the  orthogonality 
relation  [24,  Eqn.  (3.95),  p.  106], 

[R  pJQ{k'p)Jo{kp)dp  =  ~^-[Ji(kR)}26k'k, 

Jo  2* 

determine  the  constants  in  the  expression  for  <j>h , 

Ak'  +  Bv  =  2  [*  pM0)Jo(kp)dp/R2J*(kR) 

Jo 

and  k\Ak'Ck'L  -  Bk'e~k'L)  =  2  T pML)Jo(kp)dp/R2J^(kR) 

Jo 

B.4.3  Matrix  Solution 

Boundary  conditions 

We  may  specify  either  Dirichlet  or  Neumann  boundary  conditions  or  equally 
well  a  combination  of  the  two-“mixed  boundary  conditions”. 
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Dirichlet  Boundary  conditions  Here  the  value  of  the  potential  is  given 
at  a  boundary.  Assume  that  the  nth  grid  line  contains  the  given  values 
<f>i,n  =  /«•  The  differential  operator  at  the  n  —  lai  grid  line  is: 


^»',n  4~  $i,n— 2 

A2 

2<^i,n-l  "I"  2 

or - A2 - 


—  Pi,n-l 


—  Pi,n— 1 


The  boundary  value  is  absorbed  into  the  charge. 

Neumann  Boundary  Conditions  Here  the  value  of  the  electric  field 
is  given  at  a  boundary.  Assume  that  the  nth  grid  line  contains  the  given 
values  #f„  =  fif,  =  *,'ntl2~*‘‘n~1 

Then  the  differential  operator  on  the  nth  grid  line  becomes 

<t>i,n+l  ~  2 </>,>  +  <Ai',n-l  _ 

A2  _  Pi,n 


~2  <j>i<n  +  2<j)itn^l 

A2 


-  Pi,  n  A 


Again  the  boundary  value  is  absorbed  into  the  charge  so  that  we  are  free 
to  consider  homogeneous  boundary  values,  <t>h  =  0  or  (f>’h  =  0,  only. 

Mixed  Boundary  Conditions  Some  combination  of  <j>  and  <f>'  are  given 
on  separate  (therefore  not  Cauchy)  boundaries. Assume  we  have  Dirichlet 
boundary  conditions  <f>  =  0  at  z=0  and  Neumann  boundary  conditon  <f>'  =  0 
at  x=L.  The  boundary  conditions  may  remain  unspecified  in  the  r  direction. 


Difference  Equations 

To  obtain  the  discrete  analog  to  the  V2  operator  we  use  the  3-point  two 
dimensional  difference,  that  is,  since 


Di+i4>  = 

DU  = 


fii+l  ~  4>i 
A  ’ 

4>i+\  —  2<^i  +  <t>i~  i 
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The  difference  equation  is  readily  obtained  as 

1  [rm+%((t>m+l,n  ~  ^m,n)  —  ~  ^m-l,n)] 

rm  A  r2 

(^m,n+ 1  ^m,n)  (^m.n  ^m.n— l) 

A^ 


+ 


=  -47r/>m,n  (B.5) 


for  the  interior  points. 

At  the  boundaries  (particularly  at  r=0)  it  is  useful  to  apply  the  integral 
form  of  Gauss’  Law,  Eqn.  (3.1(c)). 

J  V  •  Edx  =  ^  4x  pdx 

J  E  •  dS  =  4t rq 

where  the  surface  and  volume  integrals  are  over  the  cells  of  interest. 

At  m=0 

Ar  Ar 

+  ti-ViMy)2]  =  47r9°." 

,(<fro,n  ~  <ftl,n)  [(^o.n  ~  <fto,n+l)  ~  ($0,r»-l  “  <£o,n)]  _  47T<7o,n 

*  i  «-»  T 


Ar2 


Az2 

6^0, n  —  4<£t,n  —  ^0,n+l  ~  <)> 0,n-l 


A2 


AF 

=  47r/90,n  (B.6) 


for  Ar  =  A z  =  A  and  AF  =  7r(^)2A z. 
At  m  >  0 


*[("*  +  -)2  -  (m  -  -)2]Ar2(£zm>n+i.  -  Ezm<n_t)  + 


2ir(m  +  -)ArAzi?rm+i  n  —  27r(m  —  -)ArAzf7rm_in  =  47r9m,„ 

^m,n  ^wi,n— 1  \ 

A*  ' 

r>  /  ,  *  \  *  A  ym,n  rrati,"  i  o  /  *  \  a  a  $m,n  <Am— l,n  . 

27r(m  +  -)ArAz - — - +  27r(m  - -)ArAz - — -  =  47rgm>n 


1 

2" 


2nmAr2(^ — ^±1  +  ft"’"  + 

Az  A  z  ' 

<f>m,n  ~  <t>m+l,n  ,  n_f  __  1> 

- ^ - +  2ir(m  -  - 


27rmA(2^m>n  -  <f>m,n+i  -  <t>m,n- i)  + 

27rmA[2<£m,n  -  (1  +  ^-)^m+i,»  ~  (1  -  — )^m-i,n]  =  47 Tqm,n 
Im  Z  m 

~(1  ~  2m  )<t>rn-\,n  4~  4^m>n  —  (1+  2^)<^m+l,n  ~  ~  1 


A2 


=  47rpm>n 
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Figure  B.4:  This  is  Gauss’  Law  for  m  =  0. 

for  A r  =  Az  and  AV  =  2/KmAriAz. 

For  Er  =  0  at  r  =  R  =  MAr,  we  approximate  =  ^M,n 

(3  —  ~ )<£a/,ti  -  (1  ~  ^ —)<l>M-l,n  -  <t>M, n+l  “  4>M,n-\  =  4* Pm,n 
For  <^  =  0  at  r  —  R  =  (M  -f  1)A r 


(B.7) 


4^M,n  -  (1 


2M 


)<t>M-\,n  ~  <i>M,n+\  ~  <f>M,n- 1  =  PM,n 


(B.8) 


B.4. 4  Solution  for  Electric-field 

Mixed  Boundary  Conditions  in  z 

We  use  m  =  1  at  r  =  0.  Using  E  =  —V<f>  at  interior  points  ( ngr  —  1  > 
nr  >  2,ng  —  1  >  nz  >  l, ngr  =  M  +  1): 

o  _  0m+l,n  l,7i 

•C'rm.n  —  0 


p  _  (^m,n+ 1  ^rn.n-l 

£>zm,n  —  n 
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For  ngr  —  1  >  nr  >  2,  nz  =  0 

■E’zm.O  =  Ezm, 2  —  8npm,l  +  2  X  [— - - ^(^m+l.l  ~  <^m,l) 

m  —  1 

m  -  % ,  .  , 

+  ~  T^m,  1  —  1,1 )] 

m  —  1 

=  EZm,2  —  8npmii  —  2[1  +  2^m  _  Yy]^n»+l.l 

Erm  0  —  0  at  r=R  for  no  surface  charge. 

At  r=R  <}>  =  0  for  <f>(R)  =  0 

-SzM.n  =  0 

jp  _  Tp 

Z-'zm, n+j  ^2111,11-  j 


1  (m  -f- 1  )&Erm+x  n  ETm—iin(ni  1)A 


=  47T  pm<n 


ErM+l,n  =  “ — — ~ErM-l,n  + 

M  - 1 

(  ^-  )[87TpM,n  -  2(<^Af,„  -  ^A/.n-l)  +  2(^A/,n+1  -  <^M,n)] 

=  (1  —  -j^)Er\i-l,n  + 

(1  “  ~)[87r/3W>n  +  2(<f>Mtn+i  +  ~  ±4>M,n] 

At  the  corners: 

For  m,n=l,0 

•f'ri.o  =  0  Ez\fl  =  Ez\'2  —  8irpiti  +  4(^14  —  <^2,1) 

At  (m,n)=(l,ng), 

Erl,ng  —  0  Ez\fng  =  0 

At  r=R,z— 0  Er  =  0,  Ez  =  0(no  surface  charge). 

At  r=R,z=L,  (m,n)=(ngr,ng),  Ez  =  0 

2  1 

ErM+l,N  =  (1  -  -j^)ErM-\,N  +  (1  ~  PM,N  +  4( <t>M,N-l )  “  4>M,n)\ 
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Figure  B.6:  This  shows  the  calculation  of  E  on  the  boundaries. 


B.4.5  Charge  Sharing/Force  Interpolation 


The  weighting  scheme  chosen  in  ID  may  be  interpreted  as  a  dipole  expan¬ 
sion  of  the  charge  density  about  the  midpoint  between  cells[29].  In  two 
dimensions  a  similar  approach,  but  including  a  quadrupole  term  yields  a 
scheme  known  as  area  weighting  [5,  p.  244].  This  approach  is  roughly  equiv¬ 
alent  to  the  scheme  used  here  consolidating  the  charge  as  a  ring  of  finite 
size.  In  the  limit  of  large  r,  the  two  approaches  are  exactly  equivalent. 

This  scheme  may  be  readily  derived  by  application  of  Gauss’  Law: 

J  E  •  dA  =  47 xq 

ET+(z)(2nr+Az)  —  ET-(z)(2irr~Az)  -f 
Ez+[*(rl  -  rl))  -  Ez-[*(rl  -  rl))  =  47rp[7r(r^  -  r2_)]Az 

=  6Er(z+)  +  (l-6)Er(z.) 

Er(2irrAz )  —  Er-(2nr~Az)  +  ( Ez+  —  Ez-)[n(r 2  —  rl)]  =  47r/>r[7r(r2  —  rl)]Az 


where  we  designate  + ,-  as  the  upper  and  lower  limits  of  the  particle  extent. 
rET  =  Er_r_  +  2p[ir(r2  -  r2_)]  -  AEz{^~p  +  2tt (pr  -  p)[(r2  -  rl)] 


=  Fr_r_+( 
=  Fr_r_+( 


r2  —  r2 


xr.  AEZ, 

2  )i4’r^  ^“1  + 


r2  —  r2  ET+(2nr+  Az)  —  ET-(27rr„Az) 


■)[ 


ir(rl  -rl)Az 


]  +  ■ 


=  Er-r.  +  (3 — r-T)[r+Er+  -  r_Fr_]  +  •  - 


‘  rij.  —  rl ' 


=  *  +  £r+( 


r2  —  r2 


ri  —  r2 


r  1  —  r* 


)  +  r_Fr_(-±— r)  + 


r+  —  rl 


=  r+[6Er(r+,z+)  +  (l-6)ET(r+,Z-)]( 


2  2 

r  —  r 


r+  —  rl 

+r-[*a(r-,*+)  +  (1  -  +  o(Ar) 

r+  —  rl 


(Fr(z)  «  Fr(z_)  +  ^[Fr(z+)-Fr(z-)]/Az) 
Similarly 


Fr+(27rr+Az)-f;r_(27rr_Az)+(F,-Ft_)[7r(r2 -rl)]  =  47rp*[7r(r2  -rl)]Az 
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E2  =  Ez _  +  Az[4np2  -  2 (r+ET+  -  r_£r_)/(r+  -  rl)] 

=  E„-  +  Sz(Ez+  —  Ez _) 

=  6E2+  +  (1  —  6)E2- 

=  6[E2(r+,z+)(~^-)  +  E2(r.,z+)(^^)} 

+  (!-«)[( )^.(r+ ,  *- )  +  E*(r- ,  *- )( )] 


In  two-dimensions  (r-z)  we  may  continue  to  use  the  same  charge  sharing 
and  force  interpolation  scheme  as  in  the  ID  case  for  the  z  direction.  In  r, 
however,  we  choose  to  treat  the  particle  as  a  ring  of  uniform  charge  density, 


P  = 


^+k~rli)Az 


2*r,  AzAr 
ir r'J  Az 

•+i 


< ! 
^  2 


The  charge  density  of  the  particle  remains  uniform  but  varies  with  r  as 
shown  in  Figure  B.7. 

The  force  is  similarly  calculated  from 


F,  =  2x  /  rdrpE,  s  „E„  +  where  „„ 

J  9i  Qi 

Force  weighting  is  identical  to  charge  sharing  except  as  noted  for  cal¬ 
culation  of  B.  Since  both  the  mass  and  the  charge  are  distributed  together 
and  because  of  symmetry  the  acceleration  for  each  bit  of  charge  is  the  same 
for  the  particle  treated  as  a  whole  as  for  individual  bits. 

Care  must  be  taken  in  calculation  of  the  force  due  to  the  magnetic 
field[8].  This  force  should  act  at  the  charge  center  of  the  particle  which  is 
located  at 

j  rdrp 
Jr.  i 

'--O’ 


r2 
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=  (n  +  r«  <  \ 

More  elaborate  schemes  for  charge  weighting  and  force  calculation  may 
be  used  but  one  may  often  obtain  a  more  “accurate”  (less  noisy)  simulation 
by  instead  substituting  a  greater  number  of  particles. 
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Appendix  C 

Buneman  Algorithm  Solution 
to  the  Mixed  Boundary 
Condition  Problem 


The  Buneman  algorithm  is  a  quick  way  to  solve  Poisson’s  equation  without 
the  use  of  Fourier  transforms.  In  cylindrical  coordinates  with  r  and  z  as 
variables  we  have 

1  d  d<f>  d2<j> 
rdrrdr+dzi~  ** P 
and  its  finite  difference  form, 

~  [rm+^(^m+l,n  —  ^m,n)  ~  “  ^m-l,n)] 

'  m 

■f"  (^m.n+l  2 <f>m,n  4"  ^m.n— l)  =  47 r^7m>nA  , 

where  we  have  assumed  a  uniform  grid  spacing  Ar  =  A z  —  A. 

It  is  mathematically  convenient  to  write  the  system  of  equations  as 

T<f>  „+i  +  A<f>n  +  T<t>n-i  =  — 47rTpnA2 
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where  <j>n  = 


,  T  is  the  diagonal  matrix 


and  A  is  the  symmetric  triadiagonal  matrix 


_ 3  1 

4.2 

m  —  ^  —4  m  -f  ^ 


M-|  -4 


With  the  further  transformation  T?<£n  — ►  p.  642],  T*pn A2  — >  pn, 

'  -6  \/2 

ii  -4  J75 

T~%AT~*  A=  V 


/m(m— 1) 


/m(m+l) 


the  system  may  simply  be  written 

— *  *»♦ 

^n+l  +  A<^„  +  <f>n- 1  =  — 47T  pn 

for  <i>  =  0  at  r=R  (B.8).  Assuming  the  boundary  conditions  of  symmetry 
about  z=L  and  <^(z  =  0)  =  0,  our  system  of  equations  look  like 

A<f>  i  +  <f>2  =  pi 

$j-\  +  A4>j  -f  <f>j+ 1  =  Pj  j=2,*  •  -,N-1 
— ♦ 

2<t>N-\  +  =  Pat 

and  may  conveniently  be  written  in  block  matrix  form 


I  A  I 
I  A  I 
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where  the  corner  values  depend  on  the  given  boundary  conditions.  At  the 
end  faces,  z=0  and  L,  A<f> i  +  02  =  —Anpi  and  <f>N+n  =  <f>N-n  respectively. 
Also  0o  =  0 2N  =  0. 

One  observes  that,  for  a  grid  spacing  of  L  =  N A,  the  solution  is  imme¬ 
diately  obtained  as 

<t>N  =  -A_147rpW 

By  analogy  we  can  solve  the  entire  system  for  a  smaller  grid  spacings  by 
suitably  lumping  equations  together  and  after  obtaining  the  solution  at  a 
single  grid  line  and  use  that  solution  as  a  new  boundary  condition  to  back- 
solve  our  equations  and  obtain  solutions  at  each  of  the  other  grid  lines.  Such 
a  method  is  conceivably  less  time  consuming  than  a  brute  force  method  of 
solving  the  system  such  as  Gaussian  elimination.  In  fact  the  Buneman 
algorithm  reduces  the  necessary  calculations  by  a  factor  In  N/N. 

After  multiplying  even  equations  by  -A  and  adding  adjacent  odd  equa¬ 
tions,  we  obtain  the  reduced  set 

(21  —  A2)02  +  04  —  Pi  +  j>3  —  Api 
0j-2  +  (2/  —  A2)0j  +  $j+2  =  (Pj-i  4-  Pj+i)  —  Apj 
20,-2  +  (2 1  —  A2)0/v  =  2f>N-l  —  Apu 

where  we  have  multiplied  the  last  equation  by  -A  but  the  second  to  last  by 
-21.  This  is  equivalent  to  increasing  the  system  length  by  2  and  placing  a 
mirror  charge  distribution  on  the  opposite  side. 

This  process  may  be  repeated  k  times  for  N  =  2k.  Each  time  we  shall 
have  a  new  matrix  =  2 1  —  (A^r  ^)2  in  a  system  of  equations  in  the 
same  form  but  half  the  size  of  the  original  set.  The  k  —  1'*  reduction  yields 

^(fc_1)02<k-l)  +  02*  =  Ptfk-1) 

202*_l  +  A^-1^02*  =  p[kk  ^ 

And  finally 

= #;> 

where  Aw  =  2  -  (A1*-11)2  and  p%  =  2 p£2l,  -  A(*-1)p^.~I>. 

A  convenient  (and  quick)  way  to  solve  this  equation  results  from  the 
factorization  of  A*rb  Note  A^  satisfies 

A(r)  =  2/-(A(r-1))2 
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which  is  a  polynomial  expression  in  A  and  I.  One  may  recognize  the  equation 
for  the  as  a  difference  equation  of  the  form 

Pi+ 1  =2  -  P- 


with  solution 

Pi  =  —  2cos2‘0 

We  can  factor  pk+ 1  by  expressing  it  in  terms  of  its  roots[21]: 
Pk+1  =  2-(-2cos2*0)2  =  0 


cos 


2  2k6  =  - 


1 

2 


cos  2  k6  =  ±1/^2 
2k6  =  (2n-l)7r/4 
7>*(2n  —  1) 


6  = 


2k+l2 


This  means  that  A W  may  be  written  as 


Aw  =  (-1)1+2*  f[(A  +  2Jcos0r) 


r=l 


and  the  solution  can  be  obtained  by  k  inversions  of  tridiagonal  matrices. 

If  applied  directly,  this  algorithm  is  unstable  but,  as  Buneman  discov¬ 
ered,  stability  may  be  acquired  by  more  subtle  treatment  of  the  charge 
densities  on  the  RHS.  Letting  yTj  =  A^pj]^  + 


Xj 


p(p  -  (A(r))  \p%r  +pJ+2r  -  ?jr)) 

tir  +  ~  2^r+1) 

Qj  ^  —  (xj-2r  +  Xj+2T) 

Pi0  +  (xi  ~  P(P) 


These  are  (11,  Eqns.  11.6  a,b,  11.7,  and  11.8]  and  is  the  Poisson  solver 
implemented  in  our  simulations. 
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